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4,  INTRODUCTION 

The  goal  of  the  proposed  researeh,  “Time-Resolved  and  Spectroscopic  Three-Dimensional 
Optical  Breast  Tomography,”  is  to  develop  a  safe  and  affordable  breast  cancer  detection 
modality  that  uses  noninvasive  near-infrared  (NIR)  light  for  imaging  and  diagnosing  tumors  in 
human  breast  with  high  resolution  and  specificity.  The  project  involves  developing  (a) 
experimental  approach  for  probing  the  specimen  using  NIR  light  and  collecting  two-dimensional 
(2-D)  images,  and  (b)  image  reconstruction  and  target  localization  algorithms  that  make  use  of 
the  experimental  data  for  generating  three-dimensional  (3-D)  tomographic  images. 

In  spite  of  a  very  late  start  owing  to  the  time  it  took  for  the  Research  Protocol  to  be  approved 
(USAMRMC  finalized  the  approval  on  December  21,  2004),  significant  progress  has  been  made 
during  the  reporting  period  covered  by  this  report. 

5.  BODY 

The  tasks  performed  and  the  progress  made  during  the  current  reporting  period  may  be  broadly 
grouped  as  follows: 

•  Development  of  experimental  arrangements, 

•  Development  of  an  approach  for  optical  tomographic  imaging  using  independent 
component  analysis  (OPTICA), 

•  Development  of  a  tomographic  image  reconstruction  approach  using  the  “round-trip 
matrix”  (RTM)  approach,  and 

•  Target  detection  and  localization  using  the  combined  experimental  and  theoretical 
approaches. 

We  provide  a  brief  outline  of  our  accomplishments  in  these  areas,  and  refer  to  appended 
publications  for  detailed  description  where  applicable. 

5.1.  Development  of  Experimental  Arrangements 

Three  complementary  experimental  arrangements  were  set  up  [Technical  Objective  (TO)  1, 
Task  1]  and  performance  tested.  These  include  a  time-resolved  transillumination  imaging 
arrangement,  a  spectroscopic  imaging  arrangement,  and  a  multi-source  illumination  and  multi¬ 
detector  signal  acquisition  arrangement. 

The  experimental  arrangement  for  time-resolved  imaging  makes  use  of  approximately  120-fs 
duration,  1-kHz  repetition-rate,  800-nm  light  pulses  form  a  Ti-sapphire  laser  and  amplifier 
system  for  probing  the  sample,  and  an  ultrafast  gated  intensified  camera  system  (UGICS)  for 
recording  2-D  images.  The  UGICS  is  a  compact  gated  image  intensifier  unit  fiber-optically 
coupled  to  a  charged-coupled  device  (CCD)  camera.  It  provides  an  electronic  gate  pulse  whose 
full-width-at-half-maximum  duration  can  be  adjusted  to  a  minimum  of  80-ps,  and  the  position  of 
the  time  gate  could  be  varied  over  a  20  ns  range  with  a  step  size  of  25  ps  (or  some  integral 
multiple  of  it).  The  signal  recorded  by  the  UGICS  at  a  particular  gate  position,  U  is  a  two- 
dimensional  (2-D)  image,  that  is,  a  2-D  intensity  distribution  I{x,  y,  tl)  formed  by  the  convolution 
of  the  transmitted  light  pulse  with  the  gate  pulse  centered  on  the  gate  position. 

The  NIR  spectroscopic  imaging  arrangement  uses  the  1200-1325  nm  output  of  a  mode- 
locked  Cr"^^:  forsterite  laser  and  regenerative  amplifier  system  to  illuminate  the  sample.  A 
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Fourier  space  gate  selects  out  a  fraction  of  the  less  scattered  image-bearing  photons.  A  camera 
lens  collects  and  directs  these  photons  to  the  sensing  element  of  a  320x240  pixels  InGaAs  NIR 
area  camera.  The  resulting  2-D  images  for  different  wavelengths  in  the  1225  -  1300  nm  range 
are  recorded  and  displayed  on  a  personal  computer,  and  stored  for  further  analysis  and 
processing. 

The  multi-source  illumination  and  multi-detector  signal  acquisition  arrangement  uses  both 
the  output  of  the  Ti:sapphire  laser  system  mentioned  above,  as  well  as,  the  continuous  wave 
(CW)  784  nm  radiation  from  a  diode  laser  delivered  by  a  200-|um  optical  fiber  for  sample 
illumination.  Multiple  source  illumination  is  realized  in  practice  by  step  scanning  the  sample 
along  the  horizontal  (x)  and  vertical  (y)  directions  across  the  laser  beam.  A  camera  lens  collects 
the  diffusely  transmitted  light  on  the  opposite  face  of  the  sample  and  projects  it  onto  the  sensing 
element  of  a  cooled  charged  couple  device  (CCD)  camera  (for  steady-state  measurements),  or  the 
UGICS  system  discussed  above  (for  time-resolved  measurements).  Each  illuminated  pixel  of  the 
CCD  camera  (or,  the  UGICS  system)  could  be  regarded  as  a  detector.  For  illumination  of  every 
scanned  point  on  the  input  surface  of  the  sample,  the  CCD  camera  (or,  the  UGICS  system) 
records  the  diffusely  transmitted  intensity  pattern  on  the  output  surface. 

5.2,  Optical  Tomographic  Imaging  using  Independent  Component  Analysis  (OPTICA) 

One  of  our  major  accomplishments  during  this  reporting  period  is  the  development  of  a 
novel  target  detection  and  localization  approach,  '  OPTICA  that  enables  3-D  localization  of 
target  with  high  level  of  accuracy  [TO  3,  TO  4,  Task  1,  Task  2,  Task  6].  The  approach  makes 
use  of  transmitted  light  signal  collected  by  multiple  detectors  following  multiple-source 
illumination  of  the  turbid  medium  containing  the  targets.  The  resulting  multiple  angular  views 
provide  robust  data  that  is  analyzed  using  the  independent  component  analysis^  (ICA)  of 
information  theory  to  determine  the  locations  of  targets  relative  to  the  medium  boundaries  with 
millimeter  accuracy.  OPTICA  also  provides  a  cross-sectional  image  of  the  target  within  the 
turbid  medium.  OPTICA  may  be  used  for  detection,  localization  and  imaging  of  absorptive, 
scattering  and  fluorescent  targets.  The  details  of  the  approach  were  published  and  are  attached  to 
this  document  as  Appendices  1,2,4, 5. 

5.3.  Optical  Tomography  using  ‘round-trip  matrix’ 

We  have  developed  a  new  numerical  algorithm  [TO  3,  Task  6]  for  obtaining  tomographic 
images  of  targets  in  turbid  medium,  and  tested  it  using  simulated  data."^  In  this  approach,  the 
transport  of  light  from  multiple  sources  through  the  turbid  medium  to  an  array  of  detectors  is 
represented  by  a  response  matrix  that  can  be  constructed  from  experimental  data.  The  ‘round-trip 
matrix’  (RTM)  is  constructed  by  multiplying  the  response  matrix  by  its  transpose  for  CW  (adjoint 
matrix  for  frequency-modulated  and  short-pulse)  illumination.  The  RTM  provides  a  mathematical 
description  of  transport  of  light  from  the  sources  through  the  turbid  medium  with  embedded 
objects  to  the  array  of  detectors  and  back,  and  is  similar  to  the  time-reversal  matrix  used  in  the 
general  area  of  array  processing  for  acoustic  and  radar  imaging.  The  vector  subspace  method  along 
with  Green’s  function  calculated  from  an  appropriate  forward  model  is  then  used  to  detect  the 
targets.  In  test  runs  using  simulated  data,"^  the  approach  was  able  to  locate  6  targets  at  different 
locations  within  a  turbid  medium  of  thickness  50/^  {It  =  transport  mean  free  path  of  the  medium) 
(detailed  in  Reference  4  attached  as  Appendix  3).  While  the  simulated  data  used  slab  geometry,  the 
approach  can  be  used  for  other  sample  geometries  and  to  different  forward  models. 
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5,4.  Target  Detection,  Localization  and  Imaging 

We  have  used  the  experimental  data  (obtained  using  the  experimental  arrangements  outlined 
in  Seetion  5.1)  along  with  the  OPTICA  algorithms  (Seetion  5.2)  on  different  types  of  samples  to 
obtain  target  loeation  and  cross  sectional  images  {Appendices  2,4  and  5),  and  demonstrated 
potential  of  both  OPTICA  {Appendix  1)  and  RTM  {Appendix  2)  approaches  using  simulated  data. 

The  samples  included  Intralipid- 10%  suspension  in  water  (concentration  adjusted  to  mimic 
the  scattering  coefficient,  absorption  coefficient,  and  anisotropy  factor  of  human  breast  tissue), 
breast  phantom  obtained  from  Professor  Hebden’s  Group  at  University  College  London  (UCL),^ 
a  small  piece  of  porcine  liver  embedded  in  porcine  tissue,  and  what  is  most  important  realistic 
breast  models  formed  using  ex  vivo  breast  tissues  obtained  from  National  Disease  Research 
Interchange  (NDRI)  [along  the  line  of  TO  1,  Task  2]. 

We  have  initiated  measurements  using  light  of  different  wavelengths  [TO  2,  Task  4]  to 
obtain  estimates  of  scattering  coefficient,  absorption  coefficient,  and  anisotropy  factor  of  breast 
tissue  specimens.  Details  of  the  results  are  presented  in  the  appended  publications,  and  key 
results  are  mentioned  in  the  following  ""Key  Research  Accomplishments'”  Section. 

6.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Developed  the  Optical  Tomography  using  Independent  Component  Analysis  (OPTICA) 
approach  for  detection,  3-D  localization,  and  cross  section  imaging  of  targets  embedded  in  a 
turbid  medium  {Appendices  1,  2,  4,5). 

•  Demonstrated  the  efficacy  of  OPTICA  in  detecting  and  locating  scattering,  absorptive,  and 
fluorescent  targets.  Key  results  include: 

(i)  Detection  and  localization  of  4  small  scattering  targets  in  a  breast  phantom  (UCL 
sample.  Reference  5),  the  weakest  one  having  a  scattering  coefficient  of  only  1 . 1 
times  that  of  the  intervening  media  (and  was  considered  to  be  “rather  unlikely  to 
be  found”^  within  mm  accuracy  (detailed  in  Appendix  2); 

(ii)  Detection  of  two  absorptive  targets  in  a  breast  tissue-simulating  model  medium  of 
thickness  50  times  the  transport  mean  free  path  {h)  (detailed 'm  Appendix  5); 

(iii)  Detection  of  a  4-mm  diameter  fluorescent  sphere  in  a  breast  tissue-simulating 
model  medium  of  thickness  50  times  the  transport  mean  free  path  (/<)  (detailed  in 
Appendix  4). 

•  Detection  and  localization  of  a  small  tumor  (-'  5  mm  x  5  mm  x  3  mm)  in  a  breast  model  (42 
mm  X  30  mm  x  33  mm)  assembled  using  ex  vivo  breast  tissues. 

•  Developed  a  new  “round-trip  matrix”  algorithm  for  obtaining  tomographic  images  of  targets 
in  turbid  medium,  and  tested  it  using  simulated  data.  The  approach  could  locate  all  6  targets 
when  simulated  data  was  used  (detailed  in  Appendix  2). 
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7.  REPORTABLE  OUTCOMES 

Journal  Articles 

1.  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen  and  R.  R.  Alfano,  “Three-dimensional  optieal  imaging 
of  objeets  in  a  turbid  medium  using  independent  eomponent  analysis:  theory  and 
simulation,”/.  Biomed.  Opt.  10,  051705  (2005). 

2.  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  M.  Brito,  and  R.  R.  Alfano,  “Three-dimensional 
optieal  tomographie  imaging  of  objeets  in  tissue-simulating  turbid  medium  using 
independent  eomponent  analysis,"  Appl.  Phys.  Lett.  87,  191112  (2005). 

Conference  Proceedings 

3.  W.  Cai,  M.  Alrubaiee,  S.  K.  Gayen,  M.  Xu,  and  R.  R.  Alfano,  “Three-dimensional  optieal 
tomography  of  objeets  in  turbid  media  using  the  ‘round-trip  matrix,”  in  Optieal 
Tomography  and  Speetroseopy  of  Tissue  VI,  edited  by  Britton  Chanee,  Robert  R.  Alfano, 
Bruee  J.  Tromberg,  Mamoru  Tamura,  Eva  M.  Sevick-Muraca,  Proeeedings  of  SPIE,  Vol. 
5693  (SPIE,  Bellingham,  WA,  2005),  pp.  4-9. 

4.  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Eluorescenee  optieal  tomography 
using  independent  eomponent  analysis  to  deteet  small  objeets  in  turbid  media,”  in  Optieal 
Tomography  and  Speetroseopy  of  Tissue  VI,  edited  by  Britton  Chanee,  Robert  R.  Alfano, 
Bruee  J.  Tromberg,  Mamoru  Tamura,  Eva  M.  Seviek-Muraea,  Proeeedings  of  SPIE  Vol. 
5693  (SPIE,  Bellingham,  WA,  2005),  pp.  221-224. 

5.  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Optieal  tomography  using 
independent  eomponent  analysis  to  detect  absorptive,  scattering,  or  fluorescent 
inhomogeneities  in  turbid  media,”  in  Optical  Tomography  and  Spectroscopy  of  Tissue  VI, 
edited  by  Britton  Chance,  Robert  R.  Alfano,  Bruce  J.  Tromberg,  Mamoru  Tamura,  Eva 
M.  Sevick-Muraca,  Proceedings  of  SPIE  Vol.  5693  (SPIE,  Bellingham,  WA,  2005),  pp. 
528-535. 

Presentations 

1.  M.  Xu,  M.  Alrubaiee,  W.  Cai,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Toward  optical  imaging 
of  small  tumors  in  breasts  using  cumulant  forward  model  and  Independent  Component 
Analysis,”  Paper  PI 4- 19  presented  at  the  Era  of  Hope,  Department  of  Defense  Breast 
Cancer  Research  Program  Meeting,  June  8-11,  2005,  Philadelphia,  PA. 

2.  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Optical  tomography  using 
Independent  Component  Analysis  for  breast  cancer  detection,”  Paper  P64-2  presented  at 
the  Era  of  Hope,  Department  of  Defense  Breast  Cancer  Research  Program  Meeting,  June 
8-11,  2005,  Philadelphia,  PA. 

3.  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Three-dimensional  localization  of 
scattering  targets  in  turbid  media  using  Independent  Component  Analysis,”  Paper  CThJ3 
presented  at  the  Conference  on  Lasers  and  Electro-Optics  (CLEO),  May  22-27,  2005, 
Baltimore,  Maryland. 

4.  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Eluorescenee  optical  imaging  in 
turbid  media  using  Independent  Component  Analysis,”  Paper  CEJl  presented  at  the 
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Conference  on  Lasers  and  Electro-Optics  (CLEO),  May  22-27,  2005,  Baltimore, 
Maryland. 

5.  W.  Cal,  M.  Alrubaiee,  S.  K.  Gayen,  M.  Xu,  and  R.  R.  Alfano,  “Three-dimensional  optical 
tomography  of  objects  in  turbid  media  using  the  ‘round-trip  matrix.”  Paper  5693-02 
presented  at  the  Optical  Tomography  and  Spectroscopy  of  Tissue  VII  Conference  of 
SPIE’s  BIOS  2005/  Photonics  West,  January  22-27,  2005,  San  Jose,  California. 

6.  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Eluorescence  optical  tomography 
using  independent  component  analysis  to  detect  small  objects  in  turbid  media.”  Paper 
5693-39  presented  at  the  Optical  Tomography  and  Spectroscopy  of  Tissue  VII 
Conference  of  SPIE’s  BIOS  2005/  Photonics  West,  January  22-27,  2005,  San  Jose, 
California. 

7.  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Optical  tomography  using 
independent  component  analysis  to  detect  absorptive,  scattering,  or  fluorescent 
inhomogeneities  in  turbid  media.”  Paper  5693-02  presented  at  the  Optical  Tomography 
and  Spectroscopy  of  Tissue  VII  Conference  of  SPIE’s  BiOS  2005/  Photonics  West, 
January  22-27,  2005,  San  Jose,  California. 

8.  CONCLUSION 

The  work  carried  out  during  this  reporting  period  shows  the  potential  for  detection  and  three- 
dimensional  localization  of  targets  (including  tumor)  within  a  turbid  medium  (including  a  model 
breast  formed  with  ex  vivo  tissues)  with  significant  accuracy  based  on  the  differences  in  the 
scattering,  absorption,  and  fluorescence  characteristics  of  the  target  (e.g.,  tumor)  and  the 
intervening  medium  (e.g.,  normal  breast  tissue). 

“No  What  Section" 

•  A  recent  study  involving  35,3 19  patients  underscores  the  influence  of  primary  tumor  location 
on  breast  cancer  prognosis,  and  makes  it  imperative  that  breast  cancer  detection  modalities 
obtain  three  dimensional  (3-D)  location  of  the  tumor  relative  to  the  axilla.^  The  current  work 
is  an  important  development  in  obtaining  noninvasive  3-D  location  of  a  tumor  within  the 
breast. 

•  The  applicability  of  OPTICA  for  scattering,  absorptive,  and  fluorescent  targets  makes  it 
versatile  since  all  three  phenomena  may  be  used  for  contrast  enhancement  between  the  tumor 
and  normal  breast  tissues.  Eluorescence-based  detection  may  require  use  of  contrast  agents, 
such  as,  molecular  beacons.  Use  of  far-red  and  NIR  native  fluorescence  wing^  is  also  a 
possibility. 

•  Three-dimensional  target  localization  will  enable  closer  probing  of  a  smaller  volume  around 
of  the  target  providing  more  details  since  smaller  pixel  size  could  be  used  without  increasing 
the  computation  time  (as  a  smaller  volume  will  be  probed). 

•  Eurther  development  of  the  OPTICA  approach  may  find  another  important  application  in 
lumpectomy,  where  the  surgeon  looks  for  the  margins  of  the  tumor  with  high  accuracy. 
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Abstract.  A  new  imaging  approach  for  3-D  localization  and  charac¬ 
terization  of  objects  in  a  turbid  medium  using  independent  compo¬ 
nent  analysis  (ICA)  from  information  theory  is  developed  and  demon¬ 
strated  using  simulated  data.  This  approach  uses  a  multisource  and 
multidetector  signal  acquisition  scheme.  ICA  of  the  perturbations  in 
the  spatial  intensity  distribution  measured  on  the  medium  boundary 
sorts  out  the  embedded  objects.  The  locations  and  optical  character¬ 
istics  of  the  embedded  objects  are  obtained  from  a  Green's  function 
analysis  based  on  any  appropriate  model  for  light  propagation  in  the 
background  medium.  This  approach  is  shown  to  locate  and  charac¬ 
terize  absorptive  and  scattering  inhomogeneities  within  highly  scatter¬ 
ing  medium  to  a  high  degree  of  accuracy.  In  particular,  we  show  this 
approach  can  discriminate  between  absorptive  and  scattering  inho¬ 
mogeneities,  and  can  locate  and  characterize  complex  inhomogene¬ 
ities,  which  are  both  absorptive  and  scattering.  The  influence  of  noise 
and  uncertainty  in  background  absorption  or  scattering  on  the  perfor¬ 
mance  of  this  approach  is  investigated.  ©  2005  Society  of  Photo-Optical  instru¬ 
mentation  Engineers.  [DOI:  10.1117/1.2101568] 
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1  Introduction 

Optical  probing  of  the  interior  of  multiply  scattering  colloidal 
suspensions  and  biological  materials  has  attracted  much  atten¬ 
tion  over  the  last  decade.  In  particular,  biomedical  optical  to¬ 
mography  and  spectroscopy,  which  has  the  potential  to  pro¬ 
vide  functional  information  about  brain  activities  and 
diagnostic  information  about  tumors  in  breast  and  prostate, 
are  being  actively  pursued.*”*^  Simultaneous  developments  in 
experimental  apparatus  and  techniques  for  object  interroga¬ 
tion  and  signal  acquisition,^’^'^’'*’*^  analytical  models  for  light 
propagation,  and  computer  algorithms  for  image 

reconstruction®’*  hold  promise  for  realization  of  these  poten¬ 
tials  of  optical  tomography. 

Multiple  scattering  of  light  in  thick  turbid  media  precludes 
direct  imaging  of  embedded  targets.  One  typically  uses  an 
inverse  image  reconstruction®’*  (HR)  approach  to  reconstruct  a 
map  of  the  optical  properties,  such  as  absorption  coefficient 
{/mj  and  scattering  coefficient  of  the  medium  by  match¬ 
ing  the  measured  light  intensity  distribution  on  the  boundary 
of  the  turbid  medium  to  that  calculated  by  a  forward  model 
for  the  propagation  of  light  in  that  medium.  The  commonly 
used  forward  models  include  the  radiative  transfer 
equation*^’^*  (RTE),  the  diffusion  approximation  (DA)  of  the 
RTE  (Refs.  6  and  8),  and  random  walk  of  photons. 


Address  all  correspondence  to  M.  Xu,  the  City  College  of  the  City  University  of 
New  York,  Department  of  Physics,  J41 9,  1 38th  Street  and  Convent  Avenue,  New 
York,  New  York  10031.  Tel.:  (212)  650-6865;  Fax:  (212)  650-5530;  E-mail: 
minxu@sci.ccny.cuny.edu 


The  inversion  problem  is  ill-posed  and  must  be  regularized 
to  stabilize  the  inversion  at  a  cost  of  reduced  resolution.*’^® 
Both  iterative  reconstruction  and  noniterative  linearized  inver¬ 
sion  approaches  have  been  used  to  solve  the  inversion  prob¬ 
lem  in  optical  tomography,  which  is  weakly  nonlinear,  with 
limited  success.  Reconstruction  of  images  with  adequate  spa¬ 
tial  resolution  and  accurate  localization  and  characterization 
of  the  inhomogeneities  remain  a  formidable  task.  The  time 
required  for  data  acquisition  and  image  reconstruction  is  an¬ 
other  important  consideration. 

In  this  paper,  we  present  a  novel  algorithm  based  on  the 
independent  component  analysis^^’^*  (ICA)  from  information 
theory  to  locate  absorptive  and  scattering  inhomogeneities 
embedded  in  a  thick  turbid  medium  and  demonstrate  the  effi¬ 
cacy  using  simulated  data.  ICA  has  been  successfully  applied 
in  various  other  applications  such  as  electroencephalogram 
and  nuclear  magnetic  resonance  spectroscopy.^*”*^  We  refer  to 
this  information  theory-inspired  approach  as  optical  imaging 
using  independent  component  analysis,  abbreviated  as,  OP¬ 
TICA.  The  novelty  of  OPTICA  over  other  ICA  applications  is 
that  OPTICA  associates  directly  the  independent  components 
to  the  Green’s  functions  responsible  for  light  propagation  in 
the  turbid  medium  from  the  inhomogeneities  to  the  source  and 
the  detector,  and  therefore  the  retrieved  independent  compo¬ 
nents  can  be  used  to  locate  and  characterize  the  inhomogene¬ 
ities. 

OPTICA  uses  a  multisource  illumination  and  multidetector 
signal  acquisition  scheme  providing  a  variety  of  spatial  and 


1083-3668/2005/10(51/051 705/1 2/$22.00  ©  2005  SPIE 


Journal  of  Biomedical  Optics 


051705-1 


September/October  2005  •  Vol.  10(5) 


Xu  et  al.:  Optical  imaging  of  turbid  media... 


angular  views  essential  for  three-dimensional  (3-D)  object  lo¬ 
calization.  Each  object  (or  inhomogeneity)  within  the  turbid 
medium  alters  the  propagation  of  light  through  the  medium. 
The  influence  of  an  object  on  the  spatial  distribution  of  the 
light  intensity  at  the  detector  plane  involves  propagation  of 
light  from  the  source  to  the  object,  and  from  the  object  to  the 
detector,  and  can  be  described  in  terms  of  two  Green’s 
functions  (propagators)  describing  light  propagation  from 
source  to  the  object  and  that  from  the  object  to  the  detector, 
respectively. 

The  absorptive  or  scattering  inhomogeneities  illuminated 
by  the  incident  wave  are  assumed  to  be  virtual  sources,  and 
the  perturbation  of  the  spatial  distribution  of  the  light  intensity 
on  the  medium  boundary  is  taken  to  be  some  weighted  mix¬ 
ture  of  signals  arriving  from  these  virtual  sources.  These  vir¬ 
tual  sources  are  statistically  independent  and  can  be  recovered 
by  ICA  of  the  recorded  data  set.  The  number  of  leading  inde¬ 
pendent  components  is  same  as  the  number  of  embedded  ob¬ 
jects.  The  location  and  characterization  of  inhomogeneities 
are  obtained  from  the  analysis  of  the  retrieved  virtual  sources 
using  an  appropriate  model  of  light  propagation  in  the  back¬ 
ground  medium. 

The  remainder  of  this  paper  is  organized  as  follows.  Sec¬ 
tion  2  presents  a  brief  introduction  to  ICA  and  reviews  the 
general  theoretical  framework  for  OPTICA.  Section  3  pre¬ 
sents  the  results  from  simulations  for  different  configurations. 
Implications  of  these  results  are  discussed  in  Sec.  4. 

2  Theoretical  Formalism 

2.1  ICA 

Blind  source  separation  is  a  class  of  problem  of  general  inter¬ 
est  that  consists  of  recovering  unobserved  signals  or  virtual 
sources  from  several  observed  mixtures.  Typically  the  obser¬ 
vations  are  the  output  of  a  set  of  sensors,  where  each  sensor 
receives  a  different  combination  of  the  source  signals.  Prior 
knowledge  about  the  mixture  In  such  problems  is  usually  not 
available.  The  lack  of  prior  knowledge  is  compensated  by  a 
statistically  strong  but  often  physically  plausible  assumption 
of  independence  between  the  source  signals.  Over  the  last 
decade,  ICA  has  been  proposed  as  a  solution  to  the  blind 
source  separation  problem  and  has  emerged  as  a  new  para- 
digm  in  signal  processing  and  data  analysis.  ’  ’  ' 

The  simplest  ICA  model,  an  instantaneous  linear  mixture 
model,^^  assumes  the  existence  of  n  independent  signals  Si{t) 
{i=l,2, ...  ,n)  and  the  observation  of  at  least  as  many  mix¬ 
tures  Xi{t)  (/=  1 ,2, . . .  ,m)  by  m^n  sensors,  these  mixtures 
being  linear  and  instantaneous,  i.e.,  Xi(t)  ='2f'j^^ajjS j{t)  for  each 
i  at  a  sequence  of  time  t.  In  a  matrix  notation, 

x(?)  =  As(t),  AeR'"''",  (I) 

where  A  is  the  mixing  matrix.  The  j’th  column  of  A  gives  the 
mixing  vector  for  the  y’th  virtual  source.  ICA  can  be  formu¬ 
lated  as  the  computation  of  an  nXm  separating  matrix  B 
whose  output 

y(?)  =  Bx(0  =  Cs(f),  B  e  R"''”',  C  ^  BA  e  R"""", 

(2) 

is  an  estimate  of  the  vector  s(t)  of  the  source  signals. 


The  basic  principle  of  ICA  can  be  understood  in  the  fol¬ 
lowing  way.  The  central  limit  theorem  in  probability  theory 
tells  us  that  the  distribution  of  independent  random  variables 
tends  toward  a  Gaussian  distribution  under  certain  conditions. 
Thus,  a  sum  of  multiple  independent  random  variables  usually 
has  a  distribution  that  is  closer  to  Gaussian  than  any  of  the 
original  random  variables.  In  Eq.  (2)  y,(r)=2^C,y5^(i),  as  a 
summation  of  independent  random  variables  Sjit),  is  usually 
more  Gaussian  than  Sjit),  while  y,(r)  becomes  least  Gaussian 
when  it  in  fact  equals  one  of  the  Sj{t).  This  heuristic  argument 
shows  that  ICA  can  be  intuitively  regarded  as  a  statistical 
approach  to  find  the  separating  matrix  B  such  that  y,(f)  is 
least  Gaussian.  This  can  be  achieved  by  maximizing  some 
measure  of  non-Gaussianity,  such  as  maximizing  kurtosis^^’^^ 
(the  fourth-order  cumulate),  of  y,(?). 

ICA  separates  independent  sources  from  linear  instanta¬ 
neous  or  convolutive  mixtures  of  independent  signals  without 
relying  on  any  specific  knowledge  of  the  sources  except  that 
they  are  independent.  The  sources  are  recovered  by  a  maxi¬ 
mization  of  a  measure  of  independence  (or,  a  minimization  of 
a  measure  of  dependence),  such  as  non-Gaussianity  and  mu¬ 
tual  information  between  the  reconstructed  sources. The 
recovered  virtual  sources  and  mixing  vectors  from  ICA  are 
unique  up  to  permutation  and  scaling. 

2.2  Optical  Imaging  Using  ICA 

The  classical  approach  to  propagation  of  multiply  scattered 
light  in  turbid  media,  which  assumes  that  phases  are  uncorre¬ 
lated  on  scales  larger  than  the  scattering  mean  free  path 
leads  to  the  RTE  in  which  any  interference  effects  are 
neglected.^^  The  RTE  does  not  admit  closed-form  analytical 
solutions  in  bounded  regions  and  its  numerical  solution  is 
computationally  expensive.  The  commonly  used  forward 
models  in  optical  imaging  of  highly  scattering  media  is®’*  the 
DA  to  RTE. 

The  approach  OPTICA  can  be  applied  to  different  models 
of  light  propagation  in  turbid  media,  such  as  the  diffusion 
approximation,  '  the  cumulant  approximation,  ’  ’  the  ran¬ 
dom  walk  model,*®’^'*  and  radiative  transfer^^'*'*  when  they  are 
linearized.  The  diffusion  approximation  is  valid  when  the  in¬ 
homogeneities  are  deep  within  a  highly  scattering  medium. 
We  discuss  only  the  formalism  of  OPTICA  in  the  diffusion 
approximation  in  this  paper. 

In  this  diffusion  approximation,  the  perturbation  of  the  de¬ 
tected  light  intensities  on  the  boundaries  of  the  medium,  the 
scattered  wave  field,  due  to  absorptive  and  scattering  objects 
(inhomogeneities)  can  be  written  as*’** 

<^sca(rrf,r,)  =  -  [  G(r^,r)^/r„(r)cG(r,rJdr 

5D(r)cVrG(rrf,r)  •  VrG(r,r,)dr,  (3) 

to  the  first  order  of  the  Born  approximation*®  when  illumi¬ 
nated  by  a  point  source  of  unit  power,  where  r^,  r,  and  are 
the  positions  of  the  source,  the  inhomogeneity,  and  the  detec¬ 
tor,  respectively;  Pa,ohj~ Pa  SD=Dg^-D  are  the 

differences  in  absorption  coefficient  and  diffusion  coefficient, 
respectively,  between  the  inhomogeneity  and  the  background; 
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c  is  the  speed  of  light  in  the  medium;  and  G(r,r')  is  the 
Green’s  function  describing  light  propagation  from  r'  to  r 
inside  the  background  turbid  medium  of  absorption  and  dif¬ 
fusion  coefficients  /Xg  and  D. 

Equation  (3)  is  written  in  the  frequency  domain  and  does 
not  explicitly  show  the  modulation  frequency  w  of  the  inci¬ 
dent  wave  for  clarity.  The  following  formalism  applies  to  con¬ 
tinuous  wave,  frequency-domain,  and  time-resolved  measure¬ 
ments.  The  time-domain  measurement  is  first  Fourier 
transformed  over  time  to  obtain  data  over  many  different  fre¬ 
quencies. 

The  Green’s  function  G  for  a  slab  geometry  in  DA  is  given 
by 


G(r,r')  =  G{p,z,z')  = 


GO 


1 

4i7D 


E 


exp(-  Kr^) 

rt 


exp(-  Kf^) 

r~k 


aj{rd)  =  PjG{ra,rj),  (6) 

where  aj  and  (Bj  are  scaling  constants  for  the  y’th  inhomoge¬ 
neity. 

Both  the  location  and  strength  of  the  y’th  object  can  be 
computed  by  a  simple  fitting  procedure  using  Eq.  (6).  'iVe 
adopted  a  least-squares  fitting  procedure  given  by: 

min  I  E)  [aj'ijlG)  -  G(r;,r,)]^ 

r, 

+  2  -  G{r^,r^)]4 .  (7) 

■•rf  J 

The  fitting  yields  the  location  and  the  two  scaling  constants 
aj  and  I3j  for  the  y’th  inhomogeneity  whose  absorption 
strength  is  then  given  by  qj=ajl3j. 


rt  =  [p^  +  {zT  z'  ±2kd)y'\  (4) 

for  an  incident  amplitude-modulated  wave  of  modulation  fre¬ 
quency  w,  where  A:=0,  ±  1 ,  ±2, . . .  ,p=[(x— x')^  +  (y 
is  the  distance  between  the  two  points  r={x,y,z) 
and  r'  =  {x'  ,y'  ,z')  projected  onto  the  xy  plane;  K=[{/uig 
-io)/c)/DY'^  is  chosen  to  have  a  nonnegative  real  part;  and 
the  extrapolated  boundaries  of  the  slab  are  located  at  z=0  and 
Z=d=L^+2ze,  respectively,  where  is  the  physical  thickness 
of  the  slab  and  the  extrapolation  length  Zg  should  be  deter¬ 
mined  from  the  boundary  condition  of  the  slab.^^”^^  Greens’ 
functions  in  Eq.  (3)  for  other  geometries  can  be  obtained  ei¬ 
ther  analytically  or  numerically.'**’’"** 


2.2.1  Absorptive  inhomogeneities 

We  hrst  consider  absorptive  inhomogeneities.  Under  the  as¬ 
sumption  that  absorptive  inhomogeneities  are  localized,  that 
is,  the  y’th  one  is  contained  in  volume  Vj  centered  at 
r)(l  ^7  ^7),  the  scattered  wave  field  in  Eq.  (3)  can  be  rewrit¬ 
ten  as 


-  <^sca(rd,rj  =  2  G{rj,rj)qiG{rj,r,) ,  (5) 

,/=i 

where  qj=SpbJjCj)cVj  is  the  absorption  strength  of  the  y’th 
inhomogeneity.  The  scattered  wave  is  in  a  form  of  an  instan¬ 
taneous  linear  mixture  of  Eq.  (1).  One  absorptive  inhomoge¬ 
neity  is  represented  by  one  virtual  source  qjG{rj,r^)  with  a 
mixing  vector  G(r^,ry). 

As  the  virtual  source  qjG{rj,rs)  at  the  y’th  inhomogeneity 
is  independent  of  the  virtual  sources  at  other  locations,  ICA 
can  be  used  with  the  observations  obtained  for  the  light  source 
at  n$>J  different  positions  to  separate  out  both  virtual  sources 
and  the  mixing  vectors^^’^^  rij{rj).  The  y’th  virtual 
source  and  the  y’th  mixing  vector  ay(r^)  provide  the 

scaled  projections  of  the  Green’s  function  on  the  source  and 
detector  planes,  G(r^,rj)  and  G(r^,rj),  respectively.  We  can 
write 


=  ajG{rj,r,), 


2.2.2  Scattering  inhomogeneities 

For  scattering  inhomogeneities,  a  similar  analysis  shows  the 
scattered  wave  can  be  written  as 


-  4ca(rrf,r,)  =  E  gzi^p^dWjgzi^r^s) 

,/=i 

f 

+  E  Pdj  cos  ejg^{rj,ra)q'jp,j  cos 
y=i 

j’ 

+  E  Pdj  sin  OdgxiT^p'^^ch'jPsj  sin 
1=1 

(8) 

where  =  SD{r^cV'j  is  the  diffusion  strength  of  the  y’th  scat¬ 
tering  inhomogeneity  of  volume  W  (y  =  1 ,2, . . . pjj 
=  p,j=[ix,-Xj)^+iy,-yj)^f'^;  0^ 

and  are  the  azimuthal  angles  of  r^-r^  and  r^—Vj,  respec¬ 
tively,  and  the  two  auxiliary  functions  are  given  by 


g±(r,r')^ 


4TrD 


k=-^  _ 


(xrt  +  1) 


exp(-  Arrj^) 

irt? 


-  (xrk  +  1) 


exp(-  /cr^.) 
(r-k? 


(9) 


and 


gz(Gr')  = 


4-770 


E  i  (z  -  z'  +  2kd){Krl  +  1) 


exp(-  Arrj^) 

iriy 


-{z  +  z'  -  2kd){Kr^.  +  1) 


exp(-  Kr,^) 
(r-k? 


(10) 


The  scattered  wave  from  one  scattering  inhomogeneity  is 
thus  a  mixture  of  contributions  from  (3/')  virtual  sources: 


q'jp.j  cos  q'jP.j  sin  0,g^{rj,r,), 

(11) 


with  mixing  vectors 
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Fig.  1  Light  intensity  on  one  side  of  the  slab  is  measured  when  a  point 
source  scans  on  the  other  side.  Two  inhomogeneities  are  placed  at 
(50,  60,  20)  and  (30,30,30)  mm  inside  the  slab. 


8z{rj,ra),  Pdj  cos  p^j  sin  , 

(12) 

where  1  ^7',  respectively.  Both  the  location  and  strength 

of  the  y’th  scattering  object  are  computed  by  fitting  the  re¬ 
trieved  virtual  sources  and  mixing  vectors  to  Eqs.  (11)  and 
(12)  using  a  least-squares  procedure,  respectively. 

There  are,  in  general,  three  virtual  sources  of  specific  pat¬ 
terns  (one  centrosymmetric  and  two  dumbbell-shaped)  asso¬ 
ciated  with  one  scattering  inhomogeneity,  whereas  only  one 
centrosymmetric  virtual  source  is  associated  with  one  absorp¬ 
tive  inhomogeneity.  This  difference  may  serve  as  the  basis  to 
discriminate  absorptive  and  scattering  inhomogeneities. 

The  only  assumption  made  in  OPTICA  is  that  virtual 
sources  are  mutually  independent.  The  number  of  inhomoge¬ 
neities  within  the  medium  is  determined  by  the  number  of  the 
independent  components  presented  in  the  multisource  mul¬ 
tidetector  data  set.  No  specific  light  propagation  model  is  as¬ 
sumed  in  this  step.  The  analysis  of  retrieved  independent 
components  from  ICA  then  localizes  and  characterizes  the 
absorptive  and  scattering  inhomogeneities  inside  the  turbid 
medium  using  an  appropriate  model  of  the  light  propagator. 
Extra  independent  components  may  appear,  depending  on  the 
level  of  noise  in  the  data.  These  components  can  be  discarded 
and  only  the  leading  independent  components  must  be  ana¬ 
lyzed  to  detect  and  characterize  inhomogeneities  of  interest. 


Fig.  2  Normalized  independent  spatial  intensity  distributions  on  the 
input  (or  source)  plane  (the  first  row),  the  exit  (or  detector)  plane  (the 
second  row),  and  the  least-squares  fitting  using  Eq.  (7)  (the  third  row). 
The  left  column  is  for  the  first  absorptive  inhomogeneity  at 
(50,60,20)  mm  and  the  right  column  is  for  the  second  absorptive 
inhomogeneity  at  (30,30,  30)  mm.  On  the  third  row,  the  horizontal 
profile  of  intensity  distributions  on  the  source  plane  (diamonds)  and 
on  the  detector  plane  (circles)  are  displayed,  and  solid  lines  show  the 
respective  Green's  function  fit  used  for  obtaining  locations  and 
strengths  of  objects.  The  noise  level  is  5%. 


scattered  wave  is  linear  with  respect  to  the  absorption  and 
diffusion  strengths,  we  also  set  the  strength  of  either  absorp¬ 
tion  or  diffusion  to  be  unity  in  simulations  for  convenience. 


3  Results 

Simulations  were  performed  for  a  50-mm-thick  slab,  as 
shown  schematically  in  Eig.  1.  The  absorption  and  diffusion 
coefficients  of  the  uniform  slab  are  1/300  and  D 

=  1/3  mm,  respectively,  close  to  that  of  human  breast  tissue."*^ 
The  incident  cw  beam  scans  a  set  of  21  X  21  grid  points  cov¬ 
ering  an  area  of  90  X  90  mm^.  The  spacing  between  two  con¬ 
secutive  grid  points  is  4.5  mm.  This  light  intensity  on  the 
other  side  of  the  slab  is  recorded  by  a  CCD  camera  on  42 
X  42  grid  points  covering  the  same  area. 

In  the  simulations  presented  in  the  following  subsections, 
we  fix  the  ratio  of  strength  of  absorption  to  that  of  diffusion  to 
be  0.01,  which  produce  perturbations  of  comparable  magni¬ 
tude  on  the  light  intensities  measured  on  the  detector  plane 
from  the  absorption  and  scattering  inhomogeneities.  As  the 


3.1  Absorptive  Inhomogeneities 

Two  absorptive  inhomogeneities,  each  of  a  unity  absorption 
strength,  are  placed  at  positions  (50,  60,  20)  and 

(30,30,30)  mm,  respectively.  Gaussian  noise  of  5%  was 
added  to  the  simulated  light  intensity  change  on  the  detector 
plane.  OPTICA  operates  on  a  noisy  scattering  wave 
-</>sca(r^,rj)[l +n(rrf,rj)],  where  n(r,;,r^)  is  a  Gaussian  ran¬ 
dom  variable  of  a  standard  deviation  0.05. 

ICA  of  the  perturbations  in  the  spatial  intensity  distribu¬ 
tions  provided  corresponding  independent  intensity  distribu¬ 
tions  on  the  source  and  detector  planes.  ICA-generated  inde¬ 
pendent  intensity  distributions  on  the  source  and  detector 
planes  are  shown  in  Eig.  2  for  the  two  absorptive  inhomoge¬ 
neities.  Eocations  of  the  absorptive  objects  are  obtained  from 
fitting  these  independent  component  intensity  distributions  to 
those  of  the  diffusion  approximation  in  a  slab  Eq.  (4)  by  the 
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Fig.  3  Normalized  independent  spatial  intensity  distribution  on  the  source  plane,  i.e.,  virtual  sources  (the  first  row)  and  on  the  detector  plane,  i.e., 
mixing  vectors  (the  second  row),  and  the  least-squares  fitting  (the  third  row).  The  first  column  is  for  the  first  absorptive  inhomogeneity  at 
(50,60,20)  mm  and  the  second  through  fourth  columns  are  for  the  second  scattering  inhomogeneity  at  (30,30,30)  mm  and  represent  the 
centrosymmetric  and  two  dumbbell-shaped  pairs  of  virtual  sources  and  mixing  vectors,  respectively.  The  dumbbell  comprises  one  bright  part  and 
its  antisymmetric  dark  counterpart.  In  the  third  row,  the  profile  of  intensity  distributions  on  the  source  plane  (diamonds)  and  on  the  detector  plane 
(circles)  are  displayed,  and  solid  lines  show  the  respective  Green's  function  fit  used  for  obtaining  locations  and  strengths  of  objects.  The  X 
coordinate  is  the  horizontal  coordinate.  The  p  coordinate  is  the  coordinate  on  the  symmetrical  line  passing  through  the  dumbbell  axis.  The  small 
dark  circular  region  appearing  near  the  right-upper  corner  of  the  normalized  independent  spatial  intensity  distribution  on  the  first  row  and  in  the 
fourth  column  is  an  artifact. 


least-squares  procedure  of  Eq.  (7).  The  first  object  is  found  at 
(50.0,60.0,20.0)  mm  and  the  second  one  at 
(30.0,30.0,30.1)  mm.  The  coordinates  of  both  objects  agree 
to  within  0.1  mm  of  their  known  locations.  The  strengths  of 
the  two  objects  are  o'!  =  1.00  and  172  =  0.99,  respectively,  with 
an  error  not  greater  than  1%  of  the  true  values. 

3.2  Discrimination  between  Absorptive  and  Scattering 
Inhomogeneities 

In  the  second  example,  one  absorptive  object  of  absorption 
strength  of  0.01  is  placed  at  (50,60,20)  mm  and  one  scatter¬ 
ing  object  of  diffusion  strength  of  negative  unity  (correspond¬ 
ing  to  an  increase  in  scattering  for  the  inhomogeneity)  is 
placed  at  (30,30,30)  mm,  respectively.  We  added  5%  Gauss¬ 
ian  noise  to  the  simulated  light  intensity  change  on  the  detec¬ 
tor  plane. 

Figure  3  shows  the  ICA-generated  Independent  intensity 
distributions  on  the  source  and  detector  planes  and  the  least- 
squares  fitting.  The  first  column  corresponds  to  the  absorptive 
inhomogeneity.  The  second  through  fourth  columns  corre¬ 


spond  to  the  scattering  object,  which  produces  one  pair  of 
centrosymmetric  and  two  pairs  of  dumbbell-shaped  virtual 
sources  and  mixing  vectors.  The  absorptive  inhomogeneity  is 
found  to  be  at  (50.2,60.3,20.2)  mm  with  a  strength  qi 
=0.0101.  The  scattering  object  produces  three  pairs  (one  cen¬ 
trosymmetric  and  two  dumbbell-shaped)  of  virtual  sources 
and  mixing  vectors  centering  around  the  position  (x,y) 
=  (30,30)  mm  (see  the  second  through  fourth  columns  in  Fig. 
3).  The  dumbbell-shaped  virtual  source  or  mixing  vector  com¬ 
prises  one  bright  part  and  its  antisymmetric  dark  counterpart. 
The  resolved  position  and  strength  of  the  scattering  object  are 
found  to  be  (30.0,30.0,30.0)  mm  and  q2  =  -0.99, 
(32.1 ,32.4,30.2)  mm  and  q2  =  -0.96,  and 
(31.3,30.2,27.1)  mm  and  172  =  -1.05,  respectively,  through 
fitting  to  the  Individual  pair.  For  the  scattering  object,  the  best 
result  is  obtained  from  the  fitting  to  the  first  pair  of  cen¬ 
trosymmetric  virtual  source  and  mixing  vector  from  the  scat¬ 
tering  object.  Taking  the  position  and  strength  of  the  scatter¬ 
ing  object  to  be  that  from  fitting  the  centrosymmetric  virtual 
source  and  mixing  vector,  the  error  of  the  resolved  positions 
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Fig.  4  Normalized  independent  spatial  intensity  distribution  on  the  source  plane  (the  first  row)  and  on  the  detector  plane  (the  second  row),  and  the 
least-squares  fitting  (the  third  row)  for  one  inhomogeneity  located  at  (30,30,20)  mm  with  strengths  of  absorption  =0.01  and  diffusion  92  =  1  ■  The 
first  and  second  columns  correspond  to  the  pairs  of  dumbbell-shaped  virtual  sources  and  mixing  vectors  produced  by  its  scattering  component.  The 
third  column  corresponds  to  its  absorptive  component  obtained  by  first  removing  the  scattering  contribution.  In  the  third  row,  the  profile  of  intensity 
distributions  on  the  source  plane  (diamonds)  and  on  the  detector  plane  (circles)  are  displayed,  and  solid  lines  show  the  respective  Green's  function 
fit  used  to  obtain  locations  and  strengths  of  objects.  The  X  coordinate  is  the  horizontal  coordinate  and  the  p  coordinate  is  the  coordinate  on  the 
symmetrical  line  passing  through  the  dumbbell  axis. 


of  both  objects  is  within  0.3  mm  of  their  known  locations. 
The  error  of  the  resolved  strengths  of  both  objects  is  approxi¬ 
mately  1%  of  the  tme  values. 

3.3  Colocated  Absorptive  and  Scattering 
Inhomogeneities 

For  one  complex  Inhomogeneity  that  is  both  absorptive  and 
scattering,  the  two  pairs  of  dumbbell-shaped  virtual  sources 
and  mixing  vectors  produced  by  its  scattering  abnormality  can 
be  used  to  obtain  its  scattering  strength.  By  subtracting  the 
scattering  contribution  off  the  measured  scattered  wave,  our 
procedure  can  be  applied  again  to  the  cleaned  data  and  we 
proceed  to  obtain  its  absorption  strength.  The  third  example 
considers  a  complex  inhomogeneity  at  (30,30,20)  mm  with 
strengths  of  absorption  ^i=0.01  and  diffusion  ^/2=l  (corre¬ 
sponding  to  a  decrease  in  scattering),  respectively.  We  added 
5%  Gaussian  noise  to  the  simulated  light  Intensity  change  on 
the  detector  plane. 


Figure  4  shows  the  ICA-generated  independent  intensity 
distributions  on  the  source  and  detector  planes  and  the  least- 
squares  fitting.  The  first  and  second  columns  correspond  to 
the  pairs  of  dumbbell-shaped  virtual  sources  and  mixing  vec¬ 
tors  produced  by  its  scattering  component.  The  position  and 
strength  of  this  diffusive  component  is  obtained  to  be 
(32.7,33.0,20.5)  mm  and  52=0.95,  and 

(31.7.30.1.20.4)  mm  and  52=0-96  by  fitting  the  two  indi¬ 
vidual  dumbbell-shaped  pairs,  respectively.  The  position  and 
strength  of  the  diffusive  component  is  found  to  be 

(30.9.30.9.20.4)  mm  and  52=0.95  if  both  dumbbell-shaped 
virtual  sources  and  mixing  vectors  are  used  in  fitting.  The 
third  column  of  Fig.  4  corresponds  to  its  absorptive  compo¬ 
nent  obtained  by  first  removing  the  scattering  contribution 
from  the  measured  scattered  wave.  The  depth  and  strength  of 
the  absorption  component  is  found  to  be 
(30.8,30.7,32.7)  mm  and  51=0. 0091. 
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Fig.  5  Same  as  Fig.  3  with  (a)  10%  and  (b)  20%  Gaussian  noise.  The  dumbbell-shaped  virtual  source  on  the  source  plane  in  the  fourth  column  of 
(b)  is  deformed  and  the  least-squares  fitting  is  not  shown. 
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Table  1  Comparison  of  known  and  OPTICA  determined  positions  and  strengths  of  absorption  (Abs)  and  scattering  (Sea)  inhomogeneities  under 
different  levels  of  Gaussian  noise. 


Noise 

Target 

Known  Position 
(x,y,z)  (mm) 

Known 

Strength 

Resolved  Position 
(x,y,z)  (mm) 

Resolved 

Strength 

Error  in 
Position  (mm) 

Error  in 
Strength  (%) 

5% 

Abs 

(50,  60,  20) 

0.01 

(50.2,  60.3,  20.2) 

0.0101 

-0.3 

-1 

Sea 

(30,  30,  30) 

-1 

(30.0,  30.0,  30.0) 

-0.99 

-0.1 

-1 

10% 

Abs 

(50,  60,  20) 

0.01 

(50.2,  60.3,  20.1) 

0.0101 

-0.3 

-1 

Sea 

(30,  30,  30) 

-1 

(30.0,  30.1,  30.0) 

-0.98 

-0.1 

-2 

20% 

Abs 

(50,  60,  20) 

0.01 

(50.1,  60.3,  20.5) 

0.0102 

-0.5 

-2 

Sea 

(30,  30,  30) 

-1 

(28.9,  27.0,  32.9) 

-0.59 

-3 

-50 

The  error  in  positioning  the  scattering  component  is  less 
than  1  mm  and  the  error  of  the  resolved  strength  of  the  scat¬ 
tering  strength  is  —5%.  The  errors  in  positioning  and  the 
resolved  strength  of  the  absorptive  component  equal  ~  3  mm 


and  ~10%,  respectively,  which  are  larger  than  those  of  the 
scattering  component  because  the  error  is  amplified  when  the 
estimated  scattering  component  is  used  in  subtraction  off  its 
contribution  to  the  scattered  wave  in  our  procedure. 


Fig.  6  Same  as  Fig.  3  with  (a)  40,  (b)  34,  and  (c) 


10  dB  background  absorption  uncertainty. 
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Table  2  Comparison  of  known  and  OPTICA  determined  positions  and  strengths  of  absorption  (Abs)  and  scattering  (Sea)  inhomogeneities  under 
different  levels  of  background  absorption  uncertainty. 


SNR  (dB) 

Target 

Known  Position 
(x,y,z)  (mm) 

Known 

Strength 

Resolved  Position 
(x,y,z)  (mm) 

Resolved 

Strength 

Error  in 
Position  (mm) 

Error  in 
Strength  (%) 

40 

Abs 

(50,  60,  20) 

0.01 

(50.2,  60.3,  20.1) 

0.0101 

-0.3 

-1 

Sea 

(30,  30,  30) 

-1 

(30.1,  30.1,  30.0) 

-0.99 

-0.1 

-1 

34 

Abs 

(50,  60,  20) 

0.01 

(50.2,  60.3,  20.1) 

0.0100 

-0.3 

-1 

Sea 

(30,  30,  30) 

-1 

(31.6,  31.7,  25.3) 

-0.52 

-5 

-50 

10 

Abs 

(50,  60,  20) 

0.01 

(50.6,  60.3,  20.3) 

0.0090 

-0.3 

-10 

Sea 

(30,  30,  30) 

-1 

(31.7,  29.6,  31.7) 

-0.78 

-5 

-50 

3.4  Effect  of  Noise 

To  demonstrate  the  effect  of  noise  on  the  performance  of  OP¬ 
TICA,  different  levels  of  Gaussian  noise  were  added  to  the 
simulated  light  intensity  change  on  the  detector  plane. 

Figure  5  shows  the  case  presented  In  Fig.  3  of  Sec.  3.2 
with,  instead,  10  and  20%  Gaussian  noise  added  to  the  scat¬ 
tered  wave.  The  resolved  absorptive  inhomogeneity  is  at 
(50.2,60.3,20.1)  mm  with  strength  0.0101  at  10%  noise,  and 
at  (50.1,60.3,20.5)  mm  with  strength  0.0102  at  20%  noise. 
The  resolved  position  and  strength  of  the  scattering  object  are 
found  to  be  (30.0,30.1 ,30.0)  mm  and  q2=-0.98, 

(32.1 .32.4.30.4)  mm  and  ^2=“0.95,  and 

(31.4.30.1 .27.5)  mm  and  q2=-l.00,  respectively,  through 
fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of 
dumbbell-shaped  virtual  sources  and  mixing  vectors  [see  the 
second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  10% 
noise.  The  resolved  values  become  (28.9,27.0,32.9)  mm  and 
^2  =  -0.59  from  fitting  the  pair  of  centrosymmetric  virtual 
source  and  mixing  vector  [see  the  second  column  of  Fig. 
5(b)],  and  (30.3,32.3,26.6)  mm  and  q2=-l.33  from  fitting 
the  first  pair  of  dumbbell-shaped  virtual  source  and  mixing 
vector  [see  the  third  column  of  Fig.  5(b)],  respectively,  at  20% 
noise.  The  dumbbell-shaped  virtual  source  In  the  source 
plane,  of  the  second  pair  of  dumbbell-shaped  virtual  source 
and  mixing  vector,  is  deformed  and  the  fitting  is  not  shown 
[see  the  fourth  column  of  Fig.  5(b)].  The  deformation  of 
dumbbell  appears  first  in  the  source  plane  with  the  increase  of 
noise  as  the  grid  spacing  on  the  source  plane  is  larger  than 
that  in  the  detector  plane  in  the  simulation. 

The  error  in  localization  and  characterization  of  scattering 
inhomogeneities  increases  rapidly  with  the  increase  of  noise, 
from  ~0.1  mm  in  positioning  and  ~2%  in  strength  at  10% 
noise  to  ~3  mm  in  positioning  and  ~50%  in  strength  at  20% 
noise.  On  the  other  hand,  the  effect  of  noise  on  localization 
and  characterization  of  absorptive  inhomogeneities  is  much 
smaller,  the  errors  at  both  noise  levels  are  less  than  0.5  mm  in 
positioning  and  ~2%  in  strength.  The  results  in  Sec.  3.2  and 
this  section  are  summarized  in  Table  1 . 

3.5  Effect  of  Uncertainty  in  Background 

In  the  examples  discussed  here,  we  have  assumed  the  light 
intensities  change  measured  on  the  detector  plane  is  obtained 
with  an  exact  knowledge  about  the  background.  To  examine 


the  effect  of  uncertainty  in  background  optical  property  on  the 
performance  of  OPTICA,  we  model  the  error  in  the  estimation 
of  the  background  absorption  or  diffusion  coefficients  as  a 
uniform  Gaussian  random  field  /(r) .  The  Gaussian  noise  ad¬ 
dressed  in  Sec.  3.4  is  set  to  be  zero  here.  OPTICA  operates  on 
a  “dirty”  scattering  wave  -^sca(j'd>G)  +  ^'/’sca(rd,rs),  where 
d^^s)  is  ths  change  in  the  scattered  wave  from  that  of 
a  uniform  background  of  absorption  /U,„  (or  diffusion  D)  to 
that  of  a  background  of  absorption  /x„+f{r)  [or  diffusion  D 
-l-/(r)].  The  magnitude  of  the  background  uncertainty  is  rep¬ 
resented  by  the  signal-to-noise  ratio  (SNR)  defined  by 


SNR(dB)  =  lOlogio 


E,  E,  ksca(rd,r,)p 

‘■d  *.Y _ 

^d 


(13) 


Figures  6(a)-6(c)  show  the  case  presented  in  Fig.  3  of  Sec. 
3.2  with  40,  34,  and  10  dB  SNR  due  to  background  absorp¬ 
tion  uncertainty,  respectively.  The  resolved  absorptive  inho¬ 
mogeneity  is  at  (50.2,60.3,20.1)  mm  with  strength  0.0101  at 
40  dB  SNR,  and  at  (50.2,60.3,20.1)  mm  with  strength 
0.0100  at  34  dB  SNR,  and  at  (50.6,60.3,20.3)  mm  with 
strength  0.0090  at  10  dB  SNR. 

The  resolved  position  and  strength  of  the  scattering  object 
are  found  to  be  (30.1 ,30.1 ,30.0)  mm  and  172  =  -0.99, 
(32.1 ,32.9,30.0)  mm  and  q2  =  -0.95,  and 

(31.4,30.0,27.5)  mm  and  q2  =  -l.0l,  respectively,  through 
fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of 
dumbbell-shaped  virtual  sources  and  mixing  vectors  [see  the 
second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  40  dB 
SNR.  The  resolved  values  become  (31.6,31.7,25.3)  mm  and 
q2=-0.52  and  (31.7,29.6,31.7)  mm  and  q2=-0.18  at  34  dB 
and  10  dB  SNRs,  respectively,  from  fitting  the  pair  of  cen¬ 
trosymmetric  virtual  source  and  mixing  vector  [see  the  second 
columns  of  Figs.  5(b)  and  5(c)].  The  dumbbell-shaped  virtual 
sources  and  mixing  vectors,  especially  the  dumbbell-shaped 
virtual  sources  on  the  source  plane  are  deformed  and  the  fit¬ 
ting  are  not  shown  [see  the  third  and  fourth  columns  of  Figs. 
5(b)  and  5(c)].  The  results  for  the  influence  of  background 
absorption  uncertainty  on  the  performance  are  summarized  in 
Table  2. 

Figures  7(a)  and  7(b)  show  the  case  presented  in  Fig.  3  of 
Sec.  3.2  with  34  and  10  dB  SNR  due  to  background  scattering 
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X  (mm)  X  (mm)  p  (rrm) 

Fig.  7  Same  as  Fig.  3  with  (a)  34  and  (b)  10  dB  background  scattering  uncertainty. 
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Table  3  Comparison  of  known  and  OPTICA  determined  positions  and  strengths  of  absorption  (Abs)  and  scattering  (Sea)  inhomogeneities  under 
different  levels  of  background  scattering  uncertainty. 


SNR  (dB) 

Target 

Known  Position 
(x,y,z)  (mm) 

Known 

Strength 

Resolved  Position 
(x,y,z)  (mm) 

Resolved 

Strength 

Error  in 
Position  (mm) 

Error  in 
Strength  (%) 

34 

Abs 

(50,  60,  20) 

0.01 

(50.1,  60.3,  20.1) 

0.0100 

-0.3 

-1 

Sea 

(30,  30,  30) 

-1 

(30.0,  30.1,  30.0) 

-0.99 

-0.1 

-1 

10 

Abs 

(50,  60,  20) 

0.01 

(49.9,  60.5,  20.1) 

0.0099 

-0.5 

-1 

Sea 

(30,  30,  30) 

-1 

(31.7,  31.1,  32.5) 

-0.75 

-2.5 

-25 

uncertainty.  The  resolved  absorptive  inhomogeneity  is  at 

(50.1 .60.3.20.1)  mm  with  strength  0.0100  at  34  dB  SNR 
and  at  (49.9,60.5,20.1)  mm  with  strength  0.0099  at  10  dB 
SNR. 

The  resolved  position  and  strength  of  the  scattering  object 
are  found  to  be  (30.0,30.1,30.0)  mm  and  q2=-Q.99, 
(32.2,33.0,30.0)  mm  and  q2=-0.96,  and 

(32.3.29.3.27.1)  mm  and  q2=-\.0'&,  respectively,  through 
fitting  to  the  pair  of  centrosymmetric  and  two  pairs  of 
dumbbell-shaped  virtual  sources  and  mixing  vectors  [see  the 
second  to  fourth  columns  of  Fig.  5(a)],  respectively,  at  34  dB 
SNR.  The  resolved  position  and  strength  of  the  scattering  ob¬ 
ject  are  found  to  be  (31.7,31.1,32.5)  mm  and  q2=-0.15, 
(30.9,31.4,27.5)  mm  and  q2=-\.0^,  respectively,  through 
fitting  to  the  pair  of  centrosymmetric  and  the  first  pair  of 
dumbbell-shaped  virtual  sources  and  mixing  vectors  [see  the 
second  to  fourth  columns  of  Fig.  5(b)]  at  10  dB  SNR.  The 
results  for  the  influence  of  background  scattering  uncertainty 
on  the  performance  are  summarized  in  Table  3. 

The  uncertainty  in  the  background  absorption  or  diffusion 
coefficients  affects  the  performance  of  OPTICA  in  a  similar 
fashion  as  the  noise  does  discussed  in  Sec.  3.4.  The  error  in 
localization  and  characterization  of  scattering  inhomogene¬ 
ities  increases  rapidly  while  the  error  In  localization  and  char¬ 
acterization  of  absorptive  inhomogeneities  only  increases 
mildly  with  the  increase  of  the  uncertainty  in  the  background 
optical  property.  The  uncertainty  in  background  scattering  has 
a  less  adverse  effect  on  the  performance  of  OPTICA  than  that 
In  background  absorption. 

4  Discussion 

The  simulational  study  of  OPTICA  presented  in  this  paper 
demonstrates  its  potential  in  optical  imaging  of  objects  in  tur¬ 
bid  media.  It  was  shown  to  be  able  to  locate  and  characterize 
absorptive  and  scattering  inhomogeneities  within  highly  scat¬ 
tering  medium.  In  particular,  OPTICA  can  discriminate  be¬ 
tween  absorptive  and  scattering  inhomogeneities  and  locate 
and  characterize  complex  inhomogeneltles,  which  is  both  ab¬ 
sorptive  and  scattering.  The  accuracy  of  localization  and  char¬ 
acterization  of  inhomogeneities  is  high.  In  the  cases  investi¬ 
gated  for  concentrated  inhomogeneities  within  a  tissue¬ 
emulating  slab  of  thickness  of  50  mm,  the  errors  in  resolved 
object  locations  are  not  greater  than  1  mm  and  the  errors  in 
the  resolved  optical  strengths  are  ~2%  under  favorable  noise 
levels  and  reliable  background  estimations. 


Noise  at  higher  levels  and/or  larger  uncertainty  in  the  op¬ 
tical  property  of  the  background  medium  makes  it  difficult  to 
discriminate  between  absorptive  and  scattering  inhomogene¬ 
ities.  In  such  a  situation,  other  corroborative  evidence,  such  as 
multiwavelength  measurements,  are  required  to  determine  the 
nature  of  inhomogeneities.  Noise  at  higher  levels  and/or  larger 
uncertainties  in  the  optical  property  of  the  background  me¬ 
dium  also  introduces  larger  errors  in  localization  and  charac¬ 
terization  of  scattering  inhomogeneities.  The  accuracy  of  lo¬ 
calization  and  characterization  of  absorptive  inhomogeneities 
is  only  mildly  affected  by  the  amount  of  noise  and/or  uncer¬ 
tainty  in  the  range  investigated. 

OPTICA  unmixes  inhomogeneities  based  on  the  mutual 
statistical  independence  between  them  and  does  not  rely  on  a 
Gaussian  distribution  of  light  intensity  on  the  surface  of  the 
embedding  medium.  OPTICA  has  several  salient  features. 
First,  OPTICA  provides  the  independent  components  due  to 
the  inhomogeneities  with  minimal  processing  of  the  data  and 
does  not  have  to  resort  to  any  specific  light  propagation  model 
for  obtaining  this  information.  Specific  light  propagation 
models  are  needed  only  in  the  later  stage  to  determine  loca¬ 
tion  and  optical  strength.  Second,  different  geometries,  or 
even  an  arbitrary  shaped  boundary,  can  be  used  with  OP¬ 
TICA.  Although  we  used  the  slab  geometry  in  the  work  re¬ 
ported  in  this  paper,  the  approach  does  not  depend  on  any 
specific  geometry.  Third,  the  approach  is  fast  and  is  amenable 
to  near  real-time  detection  and  localization  of  objects  in  a 
turbid  medium,  which  is  a  key  consideration  for  in  vivo  medi¬ 
cal  imaging. 

As  is  well  known,  the  diffusion  approximation  to  RTF, 
which  is  widely  used  In  Inverse  image  reconstruction,  does 
not  apply  when  the  separation  between  any  two  of  the  source, 
the  inhomogeneity  and  the  detector  is  small,  or  when  there  are 
clear  regions  in  the  medium.  A  special  treatment  is  also  re¬ 
quired  when  the  medium  has  aligned  microstructures,  such  as 
myofibrils,  axons,  or  collagen  fibers  In  tissues.^*^  The  fact  that 
a  prior  assumption  of  a  specific  model  of  light  propagation  in 
the  medium  is  not  assumed  in  the  identification  of  indepen¬ 
dent  components  by  ICA  and  is  required  only  In  a  Green’s 
function  analysis  of  the  retrieved  independent  component  is 
desirable,  especially  in  situations  that  demand  a  more  com¬ 
plex  model  than  the  conventional  DA.  Performing  the  fitting 
procedure  for  each  identified  independent  component  is  much 
simpler  and  more  transparent  than  matching  the  measured 
light  intensity  to  a  forward  model  iteratively.  The  quality  of 
reconstruction  of  OPTICA  is  expected  to  be  higher  than  the 
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conventional  approach  when  only  an  imperfect  forward  model 
is  available. 

OPTICA  is  most  suited  to  detect  small  objects.  Given  its 
ability  to  identify  low-contrast  small  objects,  the  approach  is 
expected  to  be  especially  useful  for  detection  of  tumors  at 
their  early  stages  of  development. 
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localization.  The  efficacy  of  the  approach  is  demonstrated  by  detecting  and  obtaining  location 
information  about  scattering  targets  embedded  in  human  breast  tissue-simulating  turbid  media  of 
thickness  50  times  the  transport  mean-free  path.  ©  2005  American  Institute  of  Physics. 
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Detection  and  localization  of  scattering  targets  within  a 
turbid  medium  is  a  challenging  problem  with  diverse  practi¬ 
cal  applications,  such  as  imaging  of  a  breast  tumor,  identifi¬ 
cation  of  mines  in  coastal  water,  and  detection  of  an  airplane, 
or,  structures  through  cloud  and  fog  cover.  A  recent  study 
involving  35  319  patients  underscores  the  influence  of  pri¬ 
mary  tumor  location  on  breast  cancer  prognosis,  and  makes 
it  imperative  that  breast  cancer  detection  modalities  obtain 
three  dimensional  (3D)  location  of  the  tumor  relative  to  the 
axilla.^  Optical  detection  of  targets  in  a  turbid  medium 
makes  use  of  the  difference  in  optical  properties,  such  as 
scattering  coefficient,  absorption  coefficient,  index  of  refrac¬ 
tion,  and  fluorescence  between  the  targets  of  interest  and  the 
intervening  medium.^  Multiple  scattering  of  light  by  the  tur¬ 
bid  medium  produces  a  noise  background  that  deteriorates 
the  contrast,  blurs  the  image,  and  in  severe  cases  makes  di¬ 
rect  transillumination  imaging  impossible.  Inverse  image  re¬ 
construction  approaches  that  are  commonly  used  to  retrieve 
image  information  have  to  deal  with  the  fact  that  inverse 
problems  are  ill  posed,  and  attain  different  measures  of  suc¬ 
cess  with  ~  1  cm  spatial  resolution.^ 

This  article  introduces  an  alternative  approach  for  detec¬ 
tion  and  3D  localization  of  targets  (optical  inhomogeneities) 
embedded  inside  a  highly  scattering  turbid  medium.  The  ap¬ 
proach  makes  use  of  transmitted  light  signal  collected  by 
multiple  detectors  following  multiple-source  illumination  of 
the  turbid  medium  containing  the  targets.  The  resulting  mul¬ 
tiple  angular  views  provide  robust  data  that  is  analyzed  using 
the  independent  component  analysis  (ICA)  (Ref.  4)  of  infor¬ 
mation  theory  to  determine  the  locations  of  targets  relative  to 
the  medium  boundaries  with  millimeter  accuracy.  We  refer  to 
this  optical  domain  application  of  ICA  as  OPTICA  to  distin¬ 
guish  it  from  other  applications.  While  OPTICA  may  be  used 
for  the  detection  and  localization  of  absorptive,  scattering, 
and  fluorescent  targets,  this  letter  focuses  on  scattering 
targets. 

The  perturbation  in  the  light  intensity  distribution  on  the 
boundaries  of  the  medium,  the  scattered  wave  field,  due  to 
scattering  inhomogeneities  is  given  by^ 
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^sca(i'<i,r,)  =-j  d^rSD(r)cV^G(rj,r)  ■  V,G(r,r,),  (1) 

in  the  diffusion  approximation  (DA).  Here,  r^,  r,  and  are 
the  positions  of  the  source,  the  inhomogeneity  and  the  detec¬ 
tor,  respectively,  SD=Da^-D  is  the  difference  between  the 
diffusion  coefficient  of  the  object,  Dgi,]  and  that  of  the  me¬ 
dium,  D,  c  is  the  speed  of  light  in  medium,  and  G(r ,  r^)  and 
G(r,rJ  are  source-target  and  target-detector  Green’s  func¬ 
tions,  respectively.  While  the  formalism,  detailed  elsewhere 
and  tested  for  simulated  data,®  is  applicable  for  different 
sample  sizes  and  shapes.  Green’s  functions  for  the  slab 
geometry^  are  used  here  since  rectangular  slab  samples  were 
used  in  experiments. 

Under  the  assumption  that  the  scattering  inhomogene¬ 
ities  are  localized  in  a  few  regions  within  the  turbid  medium, 
Eq.  (1)  may  be  rewritten  as 

n 

-  0sca(rrf,G)  =  E  gzi^p^dhjgzirpr,) 

Ai 

n 

+  E  Pdj  cos  9dgv{^pT^d)qjPsi  cos  6>,g^(r^-,r,) 

1=1 

n 

+  E  Pdj  sin  0dgp{rj,rd)qjPsj  sin 
1=1 

(2) 

where  qj=SD{rj)cVj  is  the  strength  of  the  jth  scattering  in¬ 
homogeneity  of  volume,  Vj  located  at  r^,  Pdj=[ixd-Xj)^ 
+  iyd-yjff'^,  Psj=\.ixs-Xjf  +  iy,-yjff'^,  and  and  0,  are 
the  polar  angles  of  r^-r^,  and  respectively,  and  gp 

are  the  longitudinal  and  transverse  components  of  the 
Green’s  functions,  and  the  summation  is  over  all  the  inho¬ 
mogeneities.  The  contribution  to  from  the  jth  inhomo¬ 
geneity  consists  of  three  terms.  These  three  terms  may  be 
interpreted  as  contributions  from  three  “virtual  sources” 
qjgzi^j’^s),  qjPsjCos  0,gp{rj,rX  and  qjp.jsin  0,gp{rj,rj) 
weighted  by  g^(ry,rj,  Pdj  cos  0dgv{rj,rf),  and 
Pdj  sin  0dg virj,rd),  respectively.  The  first  virtual  source, 
qjg^irj,rj)  is  centrosymmetric,  the  other  two  are  dumbbell 
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FIG.  1.  (a)  A  schematic  diagram  of  the  experimental  arrangement.  Inset 
shows  a  2D  array  of  horizontal  and  vertical  points  in  the  input  plane  that  ai‘e 
scanned  across  the  laser  beam,  (b)  A  typical  raw  image  recorded  by  the 
CCD.  and  how  it  is  cropped  and  binned  for  analysis. 


shaped  and  orthogonal  to  each  other.  The  centrosymmetric 
virtual  source  makes  the  dominant  contribution  to  com¬ 
pared  to  the  other  two.  The  detected  fraction  of  on  the 
exit  plane  (also  referred  to  as  the  detection  plane)  of  the 
sample  is  a  weighted  linear  mixture  of  contributions  from  3n 
virtual  sources,  if  there  are  n  scattering  targets  within  the 
medium.  All  virtual  sources  are  assumed  to  be  statistically 
independent. 

Independent  component  analysis  of  the  measured 
can  then  retrieve  the  “virtual  sources”  and  corresponding 
weighting  matrices.®  The  contribution  of  each  target  to 
on  the  detection  plane  (and  also  on  the  source  plane)  can  be 
obtained  as  a  two-dimensional  (2D)  independent  intensity 
distribution  (IID) .  IID  due  to  a  target  may  be  looked  upon  as 
the  light  intensity  pattern  that  a  source  located  at  the  target 
position  would  generate  on  the  detection  (or  source)  plane. 
The  3D  location  of  the  target  relative  to  sample  boundaries  is 
estimated  from  fits  of  these  OPTICA  generated  IID  to  the 
model  Green’s  functions. 

The  experimental  arrangement  for  demonstrating  the  ef¬ 
ficacy  of  OPTICA  is  shown  schematically  in  Fig.  1. 
Continuous- wave  784  nm  radiation  from  a  diode  laser  deliv¬ 
ered  by  a  200-/U,m-optical  fiber  was  used  for  illuminating  the 


entrance  face  (henceforth  referred  to  as  the  “source  plane”) 
of  the  slab  sample.  Multiple  source  illumination  was  realized 
in  practice  by  step  scanning  the  slab  sample  along  the  hori¬ 
zontal  (jc)  and  vertical  (y)  directions  across  the  laser  beam.  A 
camera  lens  collected  the  diffusely  transmitted  light  on  the 
opposite  face  of  the  slab  (henceforth  referred  to  as  the  “de¬ 
tection  plane”)  sample  and  projected  it  onto  the  sensing  ele¬ 
ment  of  a  cooled  charged  couple  device  (CCD)  camera.  Each 
illuminated  pixel  of  the  1024  X  1024  pixels  of  the  CCD  cam¬ 
era  could  be  regarded  as  a  detector.  For  illumination  of  every 
scanned  point  on  the  source  plane,  the  CCD  camera  recorded 
the  diffusely  transmitted  intensity  pattern  on  the  detection 
plane. 

Two  different  samples  were  used  in  the  experiments  re¬ 
ported  here.  The  first  sample  was  a  250  mm  X  250  mm 
X50  mm  transparent  plastic  cell  filled  with  Intralipid- 10% 
suspension  in  water.  The  concentration  of  Intralipid- 10%  was 
adjusted  to  provide  a  transport  length  £,  of  ~2  mm  at 
784  nm.  A  ~9-mm-diameter  glass  sphere  filled  with  a  sus¬ 
pension  of  0.707  fjLm  diameter  polystyrene  spheres  in  water 
was  the  scattering  target.  The  microspheres  do  not  absorb 
784  nm  light,  and  their  concentration  was  adjusted  to  pro¬ 
vide  a  scattering  length,  of  —0.0188  mm,  transport  length 

of  —0.133  mm,  and  anisotropy  factor,  g  — 0.858. 
The  location  of  the  center  of  the  target  was 
(25  mm, 25  mm, 21  mm)  with  respect  to  the  front  upper  left 
corner  of  the  sample  cell  [see  inset  of  Fig.  1(a)].  This  sample 
was  used  to  test  the  predictions  of  the  theoretical  formalism. 

The  second  sample  was  a  166-mm-long,  82-mm-wide, 
and  55-mm  thick  scattering  slab  cast  from  a  suspension  of 
titanium  dioxide  particles  and  a  near-infrared  dye  in  epoxy 
resin  containing  four  cylindrical  scattering  targets.*  The  slab 
material  had  an  optical  transport  length  of  —1.1  mm,  and 
absorption  coefficient  of  0.006  mm"^  Each  of  the  four  cyl¬ 
inders  had  a  length  of  5  mm,  a  diameter  of  5  mm,  absorption 
coefficient  equal  to  that  of  the  slab  material,  and  scattering 
coefficients  4,  2,  1.5,  and  1.1  times  higher  than  that  of  the 
slab.  The  center  of  each  cylinder  was  located  in  the  plane 
halfway  between  the  front  and  back  surfaces  of  the  slab,  and 
their  known  coordinates  are  presented  in  Column  3  of  Table 
I.  The  second  sample  was  used  to  test  the  efficacy  of  OP¬ 
TICA  on  a  breast-simulating  specimen.  Sample  1  was 
scanned  in  an  x-y  array  of  21  X21  grid  points  with  a  step 
size  of  2.5  mm  across  the  laser  beam,  while  Sample  2  was 
scanned  in  a  20  X  18  array  of  same  step  size. 

Figure  1(b)  presents  a  typical  2D  raw  image  of  the  de¬ 
tection  plane  recorded  by  the  CCD  camera  for  illumination 
of  a  grid  point  in  the  source  plane.  Similar  2D  raw  images 
are  recorded  for  every  scanned  grid  point  on  the  source 
plane.  As  Fig.  1(b)  further  shows,  each  raw  image  is  then 
cropped  to  select  out  the  information-rich  region,  and  binned 
to  enhance  the  signal-to-noise  ratio.  All  the  binned  images 
are  then  added  and  an  average  image  (henceforth,  referred  to 


TABLE  I.  Comparison  of  the  known  and  OPTICA  determined  positions  of  the  targets  in  Sample  2. 


Target 

No. 

/^j,target^  A^j,slab 

Known  position 

mm 

Observed  position 

mm 

Error 

(Aa',  Av,Az) 

1 

4 

(60,  60,  27.5) 

(62,  63,  28.1) 

(2,  3,  0.6) 

2 

2 

(47,  30,  27.5) 

(48,  33,  28.9) 

(1,  3,  1.4) 

3 

1.5 

(33,  60,  27.5) 

(33,  62,  27.1) 

(0,  2,  0.4) 

4 

1.1 

(20,  30,  27.5) 

(18,  33,  32.6) 

(2,  3,  5.1) 
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FIG.  2.  Independent  2D  spatial  inten¬ 
sity  distributions  at  the  detection  plane 
of  the  scattering  sphere  in  Sample  1 
generated  by  OPTICA:  (a)  The  cen- 
trosymmetric  component,  (b)  and  (c) 
the  dumbbell  shaped  components.  The 
white  dotted  circles  in  the  images  pre¬ 
sented  in  (b)  and  (c)  are  provided  as 
guide  for  the  eyes  to  show  the  high 
intensity  and  low  intensity  areas  of  the 
dumbbell.  The  white  rectangles  in  the 
images  are  the  regions  that  are  inte¬ 
grated  over  to  generate  the  spatial  pro¬ 
files.  In  the  lower  frame  of  (a),  the 
solid  line  represents  a  Green’s  func¬ 
tion  fit  to  the  experimental  data  repre¬ 
sented  by  open  squares. 


as  the  “reference  image”)  is  generated  to  serve  as  the  back¬ 
ground.  The  difference  between  the  reference  image  and  the 
mth  binned  image  is  proportional  to  the  perturbation  in  light 
intensity  distribution  on  the  detection  plane,  ^sca,m  for  illu¬ 
mination  of  the  mth  grid  point.  All  the  </>sca,m’s  are  then 
stacked,  and  used  as  input  for  independent  component 
analysis. 

OPTICA  generated  2D  independent  intensity  distribu¬ 
tions  (IIDs)  at  the  detection  plane  for  the  single  scattering 
target  in  Sample  1  are  presented  in  the  upper  frames  of  Figs. 
2(a)-2(c).  The  corresponding  spatial  intensity  profiles  inte¬ 
grated  over  the  areas  enclosed  by  white  dashed  boxes  appear 
in  the  lower  frames.  As  predicted  by  the  theoretical  formal¬ 
ism,  three  intensity  distributions  from  three  virtual  sources 
corresponding  to  the  single  scattering  target  are  observed. 
The  spatial  intensity  profile  of  Fig.  2(a)  is  symmetric  about 
the  vertical  centerline  (centrosymmetric),  the  profile  of  Fig. 
2(b)  has  a  peak  followed  by  a  dip,  and  that  of  Fig.  2(c)  has  a 
dip  followed  by  a  peak,  and  resemble  the  predicted  dumbbell 
shape.  The  location  of  the  target  determined  from  the 
Green’s  function  fit  of  the  intensity  distribution  is 
(23  mm, 25  mm, 22  mm)  which  agrees  well  with  the  known 
position  of  (25  mm, 25  mm, 21  mm).  The  relative  peak  in¬ 
tensity  of  the  centrosymmetric  component  is  approximately 
four  times  larger  than  that  of  the  dumbbells.  For  more  highly 
scattering  samples  and  with  a  decrease  in  the  signal-to-noise 
ratio,  the  dumbbell-shaped  components  get  much  reduced  in 
intensity,  and  may  not  be  observable.  It  should  be  mentioned 
that  an  absorptive  target  generates  only  a  centrosymmetric 
IID. 

The  results  of  OPTICA  measurements  on  Sample  2  are 
summarized  in  Table  I,  that  compares  the  known  locations  of 
all  4  targets  with  those  obtained  from  this  approach.  Except 
for  Target  No.  4,  coordinates  of  other  targets  are  obtained 
within  0-3  mm  (a  standard  deviation  of  ~2  mm)  of  the  re¬ 
spective  known  positions.  The  larger  error  in  the  estimated 
location  of  Target  No.  4  may  be  due  to  its  low  contrast  and 


hence  higher  susceptibility  to  noise.  Overall,  the  errors  in 
location  estimates  are  smaller  than  the  target  dimensions 
even  for  a  breast-simulating  scattering  medium. 

What  is  remarkable  about  the  OPTICA  approach  is  that 
even  at  this  initial  stage  of  development  it  could  detect  and 
locate  all  four  scattering  targets,  including  the  weakest  target 
with  a  scattering  coefficient  just  1.1  times  the  background 
and  hence  was  considered  to  be  “rather  unlikely  to  be 
found.”*  For  highly  scattering  medium  (such  as  Sample  No. 
2),  only  the  centrosymmetric  virtual  sources  corresponding 
to  scattering  targets  could  be  observed.  Additional  informa¬ 
tion,  such  as  measurements  at  different  wavelengths  are 
needed  for  identification  of  the  target  as  a  scatterer  or  an 
absorber.  Although  the  DA  was  used  in  this  case,  OPTICA 
does  not  depend  on  any  specific  light  propagation  model,  and 
can  be  used  with  other  models.  This  feature  makes  OPTICA 
a  more  general  approach  applicable  for  a  variety  of  scattering 
media  where  DA  may  not  hold.  OPTICA  is  suited  for  detec¬ 
tion  of  small  inhomogeneities  on  mm  scale,  and  has  the  po¬ 
tential  for  detection  and  localization  of  tumors  in  the  breast 
at  early  stages  of  growth. 
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ABSTRACT 

A  new  algorithm  based  on  multi-static  data  and  vector  subspace  classification  to  eigenvectors  of  a  round-trip  matrix  is 
introduced  for  optical  imaging  and  localization  of  objects  embedded  in  a  turbid  medium.  The  transport  of  light  from 
multiple  sources  through  excitation  of  the  embedded  objects  to  the  array  of  detectors  is  represented  by  a  response  matrix 
that  can  be  constructed  from  experimental  data.  The  ‘round-trip  (RT)  matrix’  is  constructed  by  multiplying  the  response 
matrix  by  its  transpose  for  continuous- wave  (adjoint  matrix  for  frequency  domain)  illumination.  Mathematically,  the  RT 
matrix  is  equivalent  to  transfer  of  light  ft'om  the  sources  via  the  embedded  objects  to  the  array  of  detectors  and  back,  and 
is  similar  to  the  time-reversal  matrix  used  in  the  general  area  of  array  processing  for  acoustic  and  radar  time-reversal 
imaging.  The  eigenvectors  with  leading  non-zero  eigenvalues  of  the  RT  matrix  correspond  to  embedded  objects,  which 
are  orthogonal  to  the  vectors  in  the  noise  subspace.  The  vector  subspace  method  along  with  Green’s  functions  calculated 
from  an  appropriate  model  for  light  propagation  through  turbid  media  is  then  used  to  determine  the  locations  of  the 
embedded  objects.  We  tested  this  algorithm  in  simulation  for  light  transmitting  through  a  50  l,r  thick  {l,r  ~  1  mm  is 
transport  mean  free  path)  parallel  slab  turbid  medium  with  up  to  six  embedded  absorptive  objects.  The  method  was  able 
to  globally  locate  all  six  objects  with  surprising  accuracy.  This  “round-trip  tomographic  imaging”  approach  is  fast, 
applicable  to  different  geometries  and  to  different  forward  models. 

Key  words:  Image  reconstruction;  Turbid  media;  Optical  medical  imaging;  Time  reversal  imaging. 

I.  Introduction 

Search  of  reliable  and  fast  running  approaches  for  inverse  image  reconstruction  of  turbid  medium  is  an  important  task  for 
optical  imaging  of  human’s  tissue,  breast,  prostate,  et  al.  Recent  inverse  algorithms,  such  as  Newton-Raphson-Marquart 
algorithms  [I]  and  that  direct  linear  inversion  of  3-D  matrices  [2],  are  time  consuming.  The  iterative  methods  may  not 
ensure  that  the  obtained  result  arrives  at  a  “global  minimum,”  and  does  not  converge  to  a  “local  minimum,”  which  is  not 
a  true  image. 

In  this  paper,  a  new  algorithm  based  on  multi-static  data  and  vector  subspace  classification  to  eigenvectors  of  a 
round-trip  matrix  is  introduced  for  optical  imaging  and  localization  of  objects  embedded  in  a  turbid  medium.  The 
transport  of  light,  from  multiple  sources  through  excitation  of  the  embedded  objects  to  the  array  of  detectors,  is 
represented  by  a  response  matrix  that  can  be  constructed  from  experimental  data.  The  ‘round-trip  (RT)  matrix’  is 
constructed  by  multiplying  the  response  matrix  by  its  transpose  for  continuous-wave  (adjoint  matrix  for  frequency 
domain)  illumination.  Mathematically,  the  RT  matrix  is  equivalent  to  transfer  of  light  from  the  sources  via  the  embedded 
objects  to  the  array  of  detectors  and  back,  and  is  similar  to  the  time-reversal  matrix  used  in  the  general  area  of  array 
processing  for  acoustic  and  radar  time-reversal  imaging.  The  time  reversal  imaging  approaches  [3-7]  were  successfully 
applied  in  areas  of  acoustic  propagation  and  coherent  propagation  of  electro-magnetic  wave  in  vacuum  space. 

The  eigenvectors  with  leading  non-zero  eigenvalues  of  the  RT  matrix  correspond  to  embedded  objects,  which 
are  orthogonal  to  the  vectors  in  the  noise  subspace.  The  vector  subspace  method,  which  is  called  “Multiple-Signal- 
Classification”  (MUSIC),  along  with  Green’s  functions  calculated  from  an  appropriate  forward  mode,  is  then  used  m 
determine  the  locations  of  the  embedded  objects.  We  test  and  apply  the  this  imaging  method  for  light  propagation  in 
highly  scattering  turbid  media,  where  the  phase  coherence  of  light  no  longer  exist,  and  photons  are  scattered  and  diffuse 
into  a  broad  area  of  detectors.  The  testing  results  using  simulated  data  show  that  this  method  is  able  to  globally  locate  Ae 
hidden  objects  with  surprising  accuracy  even  in  a  highly  scattering  turbid  medium.  When  the  locations  of  heterogeneities 
are  found,  inversion  of  optical  parameters  for  these  objects  can  be  much  easier  than  that  were  unknown.  Since  this 
imaging  method  is  simple,  running  very  fast,  and  can  be  generally  applied  for  different  geometries  and  other  conditions, 
it  is  expected  that  the  imaging  approach  have  broad  applications  in  medical  imaging 
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II.  The  steps  of  performing  imaging  using  the  round-trip  matrix 

The  following  steps  will  be  performed  to  locate  the  hidden  objects  in  a  turbid  medium, 

(1)  When  a  set  of  experimental  data  of  light  intensity  injected  from  source  Ry  and  received  by  detector  R(  ( t  7  =1, 
2,....,  N)  is  obtained,  a NxN response  matrix  K  (Kij  or  Kji )  is  constructed  by  subtracting  the  cotresponding  data  in 
a  background  medium. 

(2)  A  NxN  round-trip  matrix  T  =  K^K  for  the  frequency  domain,  or  T=  K^K  for  CW  case,  is  constructed  (either 
SDxDS  for  Ky;  or  DSxSD  for  Kji ),  where  and  are,  separately,  the  adjoint  matrix  and  the  transpose  matrix  of 
K;  The  T  matrix  is  Hermitian  (or  real  and  symmetry  in  CW  case),  which  had  N  real  eigenvalues,  and  their 
eigenvectors  are  orthogonal  with  each  other.  All  the  eigenvalues  and  the  eigenvectors  of  the  T  matrix  are 
computed  (for  example,  using  the  standard  UV  decomposition  codes). 

(3)  Separating  the  non- zero  eigenvalues  of  T  from  the  near-zero  eigenvalues,  the  number  of  non-zero  eigenvalues,  M, 
determines  the  number  of  hidden  objects.  The  corresponding  eigenvectors  are  m  =  I,....,  M  (each  has  N 
components).  Up  to  this  step,  no  physical  forward  model  of  light  propagation  is  introduced  yet. 

(4)  A  multiple  signal  classification  (MUSIC)  approach  is  used  to  determine  the  locations  of  M  hidden  objects:  The  3D 
medium  is  divided  to  p  voxels.  The  Green’s  function  is  calculated  using  a  proper  forward  model.  [G®(Xp,  Ry)  for 
SDxDS  case,  or  G^(R;,  Xp)  for  DSxSD  case],  where  Xp  is  the  position  of  pth  voxel.  These  Green’s  functions  for 
the  pth  voxel  form  a  vector  gp  of  N  components.  A  pseudo  spectrum  can  be  computed  according  to  the  following 
formula: 


S  p  ,  S  p  ^ 


M 

2 

<gp  -gp  >-  2: 

<Pm  ’gp  > 

m-l 

(1) 


where  the  inner  product  (A,  B)  =  Y.  A  (n)B(n) .  A  maximum  value  of /’(Xp)  at  Xp,  which  is  the  position  of  one 

n=l 

of  the  M  hidden  objects.  By  sorting  of  P(Xp),  the  positions  of  embedded  objects  are  determined. 


III.  Testing  based  on  simulation  data 

A  slab  turbid  medium,  with  a  transport  mean  free  path  l,r  =  1  mm,  absorption  length  4  =  300  mm,  and  thickness  Zd 
-  50  mm,  is  divided  into  20  layers.  A  CW  light  source,  injected  perpendicular  to  the  Zj  =  0  plane,  scans  by  a  16x16  array 
on  the  x-y  plane,  with  each  pixel  5x5mm^.  An  array  of  detectors  with  the  same  spacing  is  located  at  z  =  z^  plane 
(transmission  geometry).  The  medium  is  divided  into  16x16x20  voxels,  each  of  dimension  5x5x2.5  mm^.  The  indices  (i, 
j,  [)  represent  the  (x,  y,  z)  coordinates  of  the  voxel.  Six  (M  =  6)  absorbing  objects  are  located  in  the  medium,  separately, 
at  A(6,8,4),  B(13,10,10),  C(10,7,6),  D(6,9,6),  E(7,9,6),  and  F(8,12,14),  with  absorption  difference  of  App=  0.01  mm"’ 
from  the  background.  The  objects  D,  and  E  are  attached,  and  their  x-y  positions  are  close  to  that  of  A,  hence,  the  objects 
D,  E,  A,  are  difficult  to  be  distinguished.  The  objects  C,  D,  E,  are  located  in  the  same  z  layer. 

Using  a  diffusion  model  of  light  propagation  in  slab  geometry,  intensities  from  the  source  array  to  the  detector 
array  through  medium  with  and  without  hidden  objects  are  calculated  and  subtracted.  The  data  with  adding  (a)  noise  =  0, 
(b)  noise  =  5%  are  produced  as  simulated  “experimental”  data.  The  dimension  of  the  K  matrix  (sources  are  in  row,  and 
detectors  are  in  column)  is  N  =  16x16  =  256.  Then,  the  matrix  T  ~  K^K  (DSxSD)  is  constructed  from  simulated  data, 
and  256  eigenvalues  of  T  and  corresponding  eigenvectors  is  computed. 

The  computed  data  of  eigenvalues  are  sorted  in  descent  order  and  are  listed  in  the  first  column  of  the  table  1 .  At 
this  stage  (without  calculation  of  the  Green’s  function),  we  have  seen  that  there,  are  6  largest  eigenvalues,  which  can  be 
clearly  distinguished  from  other  near-zero  eigenvalues.  Hence,  the  signal  space  is  determined  as  M  =  6  from  simulated 
data.  Next  step,  the  diffusion  model  is  used  to  calculate  the  Green’s  functions  from  a  voxel  to  the  detectors  (need  not  that 
from  the  sources  to  the  voxel),  and  the  pseudo-spectrum  (16x16x20)  is  calculated  according  to  Eq.  (1)  using  these 
Green’s  function  and  the  computed  M  eigenvectors  in  the  signal  space.  The  values  of  pseudo-spectrum  are  sorted  in 
descent  order.  The  largest  values  of  them  and  the  corresponding  positions  of  voxels  are  listed  in  the  second  column  in 
the  table  1. 
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We  see  that  for  the  zero-noise  case,  the  six  non-zero  eigenvalues  are  more  10  orders  of  magnitude  larger  than  other 
near-zero  eigenvalues.  All  six  largest  pseudo-values  are  located  at  the  conect  locations  we  set  previously,  and  their 
values  are  more  10  orders  of  magnitude  larger  than  the  value  at  their  neighborhood  locations.  This  result  is  magically 

For  the  case  of  5%  noise,  the  values  of  the  near-zero  eigenvalues  greatly  increase,  but  still  one  order  of 
magnitude  less  than  smallest  non-zero  eigenvalue,  as  shown  on  the  third  column  in  the  table  1.  The  pseudo- values 
shown  in  the  fourth  column  of  the  table  1,  at  locations  (13,10,10),  (10,7,6),  and  (8,12,14),  although  have  ^eatly  reduced 
comparing  to  the  zero-noise  case,  are  clearly  maximum  comparing  to  their  neighborhood,  that  is  the  well  resolved  case 
we  discuss  before.  The  pseudo-value  at  (6,8,4)  is  also  a  maximum.  The  pseudo-values  at  (6,9,6)  and  a,9,6)  can  be 
regarded  as  a  local  maximum,  although  the  maximum  for  object  at  (7,9,6)  has  been  little  bit  moved.  All  these  peaks  are 
among  the  largest  values  in  256x20  values  in  the  pseudo-spectrum.  This  image  result  is  quite  well  companng  to  that  ot 
other  image  reconstruction  approaches. 


(a)  simulated  data  (no  noise) 


(b)  simulated  data  (5  %  noise) 


eigen  value 

S.S78267e-10 
4.244296e-ll 
1.584357e-ll 
3.861S00e-12 
1.0S489Se-13 
8.7S2975e-I5 
3.2394830-25 
7.2262780-26 
6.9364220-26 
6.6890550-26 
6.4117250-26 
6.3266710-26 
6 . 0552720-26 


i  j  I 

13  10  10 

7  9  6 

10  7  6 

6  9  6 

8  12  14 

6  8  4 


Pii  J,  0 

3.0714040-1-15 
2 .982245e-i-15 
1.27503Se-i-lS 
7 .403212e-i-14 
5.251203e-4-14 
4 .0445176-1-14 
5912 .727 
4089 . 3306 
3767 . 8683 
3053.547 
2816 .3402 
2485.2019 
2065 . 8781 


eigen  value 

5.878410e-10 
4.2394790-11 
1.S824O20-I1 
3 .8582760-12 
1.1365580-13 
-1.4730450-14 
-9.8154380-15 
-8.4442240-15 
7.1921330-15 
-6 . 4703220-15 
6.1564190-15 
-4.6794830-15 
-4.5086370-15 
4.0285350-15 
3.6575690-15 
3 .3722980-15 
3 .0879450-15 
-2.8386560-15 


i  j  I 

13  10  10 
10  7  6 

8  12  14 


13  10  9 

10  7  4 

6  9  6 

7  9  3 

6  9  7 

13  10  11 


51384.399 
23437.917 
18144.901 
2200 . 9195 
2105. 3805 
1697 . 181 
1476.2329 
1327 . 5666 
1271.6959 
977.96331 
924.28242 
923 . 62776 
791,29503 
711 . 70233 
656 . 91035 
625.15459 
596.06539 
554.31961 


IV.  Advantages  of  the  approach 

( 1 )  This  approach  is  simple  and  running  very  fast.  If  the  weight  matrix  in  the  above  example  is  written,  it  is  256x256  in 
row  and  256x20  in  column.  In  “round-trip  matrix”  imaging,  the  problem  is  reduced  to  make  a  singular  value 
decomposition  of  a  matrix  of  256  in  row  and  256  in  column.  Computation  can  be  completed  in  a  few  minute  on  a 

personal  computer.  ■*  n 

(2)  As  shown  above,  the  positions  of  objects  can  be  globally  determined,  and  the  image  results  are  quite  ' 

(3)  Only  one  kind  Green’s  function  is  needed,  either  that  from  sources  to  a  voxel  or  that  from  a  voxel  to  detectors. 
Hence,  the  requirement  for  the  forward  model  is  reduced.  Especially,  error  due  to  use  of  the  diffusion  model  in 
computation  of  the  Green’s  function  from  the  sources  to  voxels  can  be  avoided,  if  only  the  Green  s  function  trom 

voxels  to  detectors  is  required.  .. 

(4)  This  approach  can  be  generally  used  under  many  different  conditions,  either  parallel  geometry  or  cy  in 

geometry. 
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(5)  The  background  medium  may  be  uniform  or  non-uniform,  as  long  as  the  Greens  function  of  the  background  can  be 
estimated.  For  example,  this  algorithm  may  be  used  for  imaging  based  on  the  difference  of  data  obtained  using 
different  wavelengths. 


V.  Brief  description  of  theory  for  the  imaging  using  the  round-trip  matrix 

The  theoretical  foundation  of  the  above  approach  is  based  on  the  work  by  A.  J.  Devaney  [3]. 

Experimental  data  are  collected  by  a  array  of  multiple-sources  and  multiple-detectors,  which  are  expressed  by  a 
response  matrix  K.  In  a  linear  form  using  the  Bom  perturbation  method,  the  K  matrix  can  be  theoretically  expressed  as: 


K 


=  I  g‘ 

m=l 


,R 


M  ^ 
0=  ^  Sr. 

m-l 


(2) 


or 


K  =  lKj^l] 


V  ^ 

2.  8m^n 


d 

Sm 


M 
=  E 

m=l  ■'  m=l 

where  I,  j  =  1,2, . N  are,  separately,  indices  of  detectors  and  sources,  m=  1,2, . A/  is  the  index  of  hidden  objects, 

with  M  <  N;  R;,  R/,  and  X„  are,  respectively,  the  positions  of  a  source,  a  detector,  and  an  object;  r„  is  difference  of  the 
optical  parameters  (absorption  and  scattering)  of  the  object  from  that  of  the  background  medium  (the  background 
medium  may  be  uniform  or  non-uniform),  and  G^(X„„  R^)  and  G'^fR;,  X„)  are  the  Green’s  functions  in  the  background 
medium.  The  vector  =  [G(X^,Rl),G(X^,R2),-  ■  •,G(X^,R^)]  in  Eq.  (1)  has  components. 

When  data  of  the  K  matrix  is  obtained  by  measurements,  a  NxN  round-trip  matrix  T  =  K^K  for  the  frequency 
domain,  or  T=  K^K  for  CW  case,  is  constructed  (SDxDS  for  Kif,  or  DSxSD  for  Kji ),  where  and  are,  separately,  the 
adjoint  matrix  and  the  transpose  matrix  of  K.  This  T  matrix  was  called  the  “time-reversal”  matrix,  as  time-reversal 
process  in  time-domain  becomes  phase  conjugation  in  the  frequency-domain  through  a  Fourier  transform.  In  the  time- 
reversal  experiments,  a  special  “phase  conjugate”  mirror  is  applied,  that  the  process  after  reflected  from  the  mirror  is 
equivalent  to  a  time-reversal  process,  instead  of  a  process  along  future  of  time.  However,  we  do  not  see  this  time- 
reversal  process,  especially,  in  CW  case.  What  we  see  is  a  round-trip  from  sources  to  detectors,  and  then  back  to  sources 
(or,  theoretically,  from  detectors  to  source,  and  back  to  detectors).  Therefore,  a  “round-trip”  matrix  may  be  proper  for 
description  of  the  T  matrix  here.  The  T  matrix  is  Hermitian  (or  real  and  symmetry  in  CW  case),  hence,  its  eigenvalues 
are  real  and  the  corresponding  eigenvectors  consist  a  complete  set  of  orthonormal  vectors. 

In  the  many  cases,  it  can  be  proved  that  the  rank  of  T  matrix  equals  to  the  number  of  hidden  objects  M,  which 
means  this  matrix  has  M  non- zero  eigenvalues.  The  corresponding  eigenvectors  construct  a  signal  subspace.  The  N  —M 
eigenvectors  with  zero  eigenvalues  (or  near  zero,  in  the  case  of  presence  of  noise)  construct  a  noise  subspace.  These  two 
subspaces  are  orthogonal  with  each  other.  Use  of  orthogonality  of  eigenvectors  of  the  T  matrix  is  the  main  tool  for 
accurately  locating  hidden  objects  in  a  turbid  medium  using  time-reversal  imaging. 

The  T  matrix  (for  DSxSD)  can  be  written  as; 

*r- 


r  = 


M 
I  i 
Lm-\ 
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iSm 
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^m,m'Sm  Sm' 


s 
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N  * 
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,Ry)G^X, 
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(4) 


where  the  inner  product  is  defined  as  (A,  =  ZA  (n)S(n) . 

n=l 


For  a  well-resolved  case,  we  have  <  8,^,8 m 

N 

-ffi)  =  Z 


_  d  d 
'  ^  Sm'  Sm' 


n=l 


C^(R„r)*G^R„X^) 


>«  0  when  m^m' .  We  interpret  the  quantity 

(5) 


as  being  the  point  spread  function  (PSF)  of  the  detector  array,  which  represents  an  “image”  of  object  m  from 
measurements  on  N  detectors  in  terms  of  the  backprojection  of  the  outgoing  wave  Green’s  function.  If  the  object  m’  is 
located  far  enough  from  the  object  m,  that  is  X^) «  0,  then  these  two  objects  are  regarded  to  be  well  resolved. 
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If  all  M  objects  are  well  resolved  each  other,  we  have  K 


s  _  c 
mm'  ~ 


M  s  d*  d 

■  2  gffi  gfn* 

m-l 


Kn  =1  ^  Pm  >  with  a  normalized  factor  =<  Sm  >  •  I"  we  can  prove  that  M  non-zero  eigenvectors 

are  precisely  =  g£  -  and  the  corresponding  non-zero  eigensvalues  are  |  P  pj pi • 


s  d  d 

Pm^  P  m^im^ 


j*  M  \  A  A  ^  ^  s  d  d  (y  i  |2  s  d  d  /£\ 

[prove]  Tgi  =  S  A;„  ^  ^/n  Pm^m.m^  "1  1  P  m^^m^  •  ^  ^ 

The  non-zero  eigenvalue  of  mth  object  is  related  to  the  optical  strength  of  the  object,  and  the  Green’s  function  from  this 
object  to  the  sources  and  detectors.  In  the  well  resolved  case,  when  the  eigenvectors  m  the  signal  subspace  are 
computed  using  a  Green’s  function  based  a  forward  model.  The  locations  of  targets  can  be  determined  by: 


(jr^(r)=  I  G{r,^j)Pm(j)  -  ^2 

;=i  (Pm) 


H  (r,X™)  , 


where  H(r,  XJ  is  the  point  spread  function  of  a  target  located  at  X,„.  The  maximum  of  H(r,  X„)  at  r  determines  the 

positiOT  of^ith  t^gettors  subspace  Ps-  Because  T  is  a  AfxAi  matrix,  other  W-M  eigenvectors  have  zero 

eigenvalues  and  span  the  noise  subspace  Pn,  and  Ps  +  Pn  =  In  the  case  of  non-\vell-reso  ve  targe  s  e  ve 
characterizing  an  object  can  be  a  Unear  combination  of  several  eigenvectors  in  the  signal  subspace  hence  Eqjp  is  not 
available  for  determining  locations  of  objects.  A  MUSIC  approach  is  applied  to  determine  the  locations  of  winch 

uses  the  property  that  any  linear  combination  of  eigenvectors  in  the  signal  subspace  must  be  orthogonal  to  the  vectors  in 
the  noi.sc  subspace.  We  use  a  forward  model  to  calculate  the  Green’s  function  at  test  locaUon  Xpi 

g(Xp)  =  [G(Xp,Ri),G(Xp,R2),-  ■  -.GCXp.Rjv)]  .  (8) 

MUSIC  uses  a  pseudo-spectrum  to  determine  the  target  locations.  After  the  eigenvectors  in  the  noise  space  Uj  are 
obtained,  the  following  pseudo-spectrum  can  be  computed: 

P{\p)  = - ^ ^ 

SA^.«Ol<«y>g(Xp))l 

that  P(X  )  become  large  at  object  locations  and  will  be  small  otherwise.  It  is  more  convenient  to  use  P^  =I-Ps,  and 
apply  the  eigenvectors  in  the  signal  space  instead  of  the  eigenvectors  in  the  noise  space.  We  have 

z  (pj  >Sp)  =\U-P^]gp\^=  8p-  '^Pm{Mm  ^Sp)  =\8  p\  -  ^\<  Pm  ’  8  p  •  O®) 

j=M+l  \  ^  I  in=l  '  ' 

'2> 

The  denominator  of  the  Eq.  (1)  is  obtained  from  Eq.  (10).  We  use  \g p\  in  the  numerator  of  Eq.  (1)  to  make  the 
expression  normalized. 


VI.  Summary 

In  this  paper  we  test  the  image  algorithm  using  the  round-trip  matrix  to  locate  hidden  objects  i"  highly  f  4 

media  using  light  signals  of  an  array  of  multiple  sources-detectors.  The  testing  results  using  simulated  dahi  show  that  h  s 
method  is  able  to  globally  locate  the  hidden  objects  with  high  accuracy  m  a  highly  scattering  turbid  medium.  Since 
image  method  is  simple,  running  very  fast,  and  can  be  generally  applied  for  different  geometries  and  other  conditions. 

is  expected  that  the  approach  has  broad  applications  in  medical  imaging.  ,-,y-,Tvrv  I'm  a  c  a  r  rani  Mo- 
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Fluorescence  optical  tomography  using  independent  component 
analysis  to  detect  small  targets  in  turbid  media 
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ABSTRACT 

A  new  approach  for  optical  fluorescence  tomographic  imaging  of  targets  in  a  turbid  medium  that  uses  the  independent 
component  analysis  (ICA)  from  information  theory  is  presented.  Fluorescence  signals  from  targets  embedded  in  a  turbid 
medium  are  measured  on  the  boundary  of  the  medium  using  a  multi-source  excitation  and  a  multi-detector  acquisition 
scheme.  Difrercnces  between  excitation  and  fluorescence  wavelengths  enable  sensitive,  minimal-background  signal  ac¬ 
quisition.  ICA  of  the  fluorescence  signal  on  the  medium  boundary  sorts  out  the  embedded  objects,  and  their  locations  are 
obtained  from  Green’s  fiinction  analysis  based  on  any  appropriate  light  propagation  model.  Fluorescence  tomographic 
imaging  experiments  were  carried  out  using  Intralipid-1 0%  suspension  in  water  contained  in  a  50-mm  thick  rectangular 
transparent  plastic  cell  as  the  turbid  medium,  and  small  glass  spheres  containing  indocyanine  green  (ICG)  solution  as 
fluorescent  targets.  The  near-infrared  (NIR)  fluorescence  was  excited  using  785  nm  light,  and  monitored  over  a  narrow 
band  around  830  nm.  The  transport  mean  free  paths  at  785  nm  and  830  nm  were  1.01  mm  and  1.14  mm,  respectively. 
The  approach  could  image  and  determine  the  position  of  an  ICG  filled  sphere  of  radius  as  small  as  4  mm.  It  is  applicable 
to  small  objects,  different  medium  geometries,  and  amenable  to  near  real  time  imaging  applications. 

Keyword:  Near-infrared  fluorescence  imaging,  imaging  through  turbid  media,  tissue-like  phantom,  inverse  image  recon¬ 
struction,  independent  component  analysis. 

1.  INTRODUCTION 

The  advent  of  novel  fluorescence  beacons  and  contrast  agents  with  ability  to  attach  themselves  to  a  desired  large  cluster 
of  abnormal  cells  or  organelles,'  and  fluoresce  when  excited  with  light  of  appropriate  wavelength  make  fluorescence 
tomography  attractive  for  detection  and  diagnosis  of  tumors  in  human  body  organs,  such  as,  breast  and  prostate.  Direct 
imaging  approaches  fail  to  provide  accurate  location  of  fluorescent  targets  embedded  in  a  turbid  medium^’^  Researchers 
then  resort  to  tomographic  image  reconstruction  (TIR)  schemes  that  use  fluorescence  signals  measured  from  the  surface 
of  the  specimen,  an  appropriate  model  for  light  propagation,  and  numerical  algorithms  to  reconstmet  fluorescent  tar- 
get{s)  inside  the  specimen."'''' 

We  introduce  a  novel  optical  tomography  approach  that  uses  independent  components  analysis'®’”  of  information 
theory  to  analyze  the  fluorescence  signal  excited  by  a  multi-source  illumination  scheme  and  measured  by  a  multi¬ 
detector  signal  acquisition  arrangement  providing  a  variety  of  spatial  and  angular  views  essential  for  three-dimensional 
(3D)  localization  and  characterization  of  the  target.  We  refer  to  this  approach  as  fluorescence  OPTICA. 

2.  THEORY 

Fluorescence  tomography  requires  accounting  for  propagation  of  excitation  light  beams  that  induce  fluorescence  from 
embedded  targets,  and  that  of  the  fluorescence  emitted  by  the  localized  targets  to  the  detectors.  In  fluorescence  OPTICA, 
the  incident  light  at  wavelength  propagates  from  the  source  point  Ts  at  the  input  (or,  source)  plane  to  the  target  located 
at  position  r  inside  the  medium.  Fluorescence  at  wavelength  propagates  from  target  at  r  to  the  detector  at  ra.  The 
propagation  of  excitation  and  fluorescent  light  beams  is  described  by  coupled  difiusion  equations  at  the  excitation  and 
emission  wavelengths 

The  fluorescence  signal  I7„(rd,r„(o)  can  be  expressed  in  terms  of  the  two  Green’s  functions  G;c(r,rs,co)  and 
Gm(rd,r,0))  describing  light  propagation  from  the  source  at  r,  to  the  target  at  r  at  the  excitation  wavelength  and  the 
light  propagation  from  the  target  to  the  detector  at  rj  at  the  emission  wavelength  The  subscripts  x  and  m  refer  to  the 
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quantities  associated  with  the  excitation  and  emission  wavelengths,  respectively.  Assuming  the  f''  fluorescent  target  is 
contained  in  volume  Fj  centered  at  rj,  the  fluorescence  signal  under  illumination  by  a  unit  point  source  at  r,  is  given  by 

2G„(rj,ri,(0)//(o)G,(rpr„a))  (1) 

where  ~  7(ry)cf^/[l-iOJT(rj)]  represents  the  fluorescence  strength  of  the  7*  inhomogeneity,  j  is  the  fluorescent  yield, 
and  X  is  the  fluorescence  lifetime.  The  contribution  of  j*  source,  r*,  to)  is  weighted  by  a  mixing  vector 

Gm{rd,rj,a))  and  all  the  weighted  sources  are  mixed  to  produce  the  detected  signal  C4(rd,rs,a)).  ICA  assumes  that  the  sig¬ 
nals  from  the  fluorescent  targets  are  statistically  independent,  and  that  the  measured  signal  is  a  weighted  mixture  of  con¬ 
tributions  from  these  independent  sources."’’’'  By  seeking  the  maximal  mutual  independence,  both  fluorescent  sources 
and  mixing  vectors  can  be  obtained  by  independent  component  analysis  of  the  multi-source  multi-detector  data  set 

fsj  ®)* 

3.  EXPERIMENTAL  ARRANGEMENT 


Figure  1  shows  schematically  the  experimental  setup  used  for  testing  the  fluorescence  OPTICA  approach.  The  turbid 
used  in  the  experiment  was  Intralipid- 10%  suspension  in  water  contained  in  a  250  mm  x  250  mm  x  50  mm  transparent 
plastic  cell.  The  concentration  of  fntralipid-10%  was  adjusted  to  provide  a  transport  length  /<  ~lmm  at  785  ran,  The 
fluorescent  target  consists  of  glass  sphere  with  a  hook  attached  to  a  white  thin  string  as  shown  iii  Fig.  2(a).  The  sphere  is 
filled  with  a  solution  of  Indocyanine  green  (ICG)  dye  in  water.  The  absorption  coefficient  of  the  ICG  solution  at  785  nm 

is  11.5  mm'’.  To  test  how 
small  an  object  could  be 
observed,  we  used  three 
different  spheres  with  di¬ 
ameters  12  mm,  9  mm  and  4 
mm  in  turn.  A  785-mn  con¬ 
tinuous  wave  (CW)  laser 
beam  was  delivered  through 
a  ~200-|rm  diameter  optical 
fiber  and  focused  to  a  ~1 
mm  spot  onto  the  input  sur¬ 
face  of  the  sample  cell.  A 
computer-controlled  transla¬ 
tion  stage  scanned  the  sam¬ 
ple  in  a  set  of  10  x  10  grid  points  (step  size  =  2.5-mm),  in  the  x-y  plane  perpendicular  to  the  incident  beam  realizing  a 
multi-source  illumination  scheme.  For  illumination  at  everyone  of  these  points  a  1024  x  1024  pixels  cooled  charge  cou¬ 
ple  device  (CCD)  camera  acquired  a  2-D  fluorescence  image  of  the  exit  surface  (opposite  to  the  input  surface)  of  the 
cell.  The  fluorescence  spectrum  of  ICG  in  water,  when  excited  at  785-nm,  spans  the  wavelength  range  of  790-966  nm  as 
shown  in  the  inset  of  Fig.  1 .  We  used  a  narrow-band  filter  centered  at  830-nm  to  acquire  a  fraction  of  the  fluorescence, 
and  discriminate  against  the  much  more  intense  laser  excitation  beam.  The  images  were  stored  in  a  personal  computer 
for  subsequent  analysis  by  the  numerical  algorithm  of  fluorescence  OPTICA. 


Fig.l  A  schematic  diagram  of  the  experimental  arrangement  for  fluorescence  OPTICA.  The 
inset  shows  the  fluorescence  spectrum  of  ICG  solution  in  water 


4.  RESULTS 


A  typical  2-D  fluorescence  image  of  the  exit  surface  of  the  turbid  medium  recorded  by  the  CCD  camera  is  shown  in  Fig¬ 
ure  2(b).  The  fluorescent  target  in  this  case  was  the  4-mm  dicuneter  sphere.  Independent  components  analysis  used  the 
set  of  100  images  to  generate  independent  intensity  distribution  on  source  and  detector  planes,  as  shown  in  Figures  2(c) 
and  2(d),  respectively.  The  Green’s  function  (solid  line)  fit  to  these  intensity  distributions  represented  by  squares  for  the 
source  plane  and  by  circles  for  the  detector  plane  are  shown  in  Figures  3(a)  and  3(b),  respectively.  The  x-y-z  co-ordinates 
of  the  fluorescence  target  were  estimated  to  be  approximately,  x  =  15  mm,  y  =  15  mm,  and  z  =  16  mm,  which  are  in  ex¬ 
cellent  agreement  with  the  known  location  of  x  =  15  mm,  y  =  15  mm,  and  z  =  15  mm.  This  result  is  for  the  smallest  of 
the  spheres,  and  hence  the  most  difficult  of  the  three  cases.  The  locations  of  the  other  two  spheres  were  readily  observed 
using  the  fluorescence  OPTICA. 


222  Proc.  of  SPIE  Vot.  5693 


Intensity  (arb.  units) 


(c)  .  (<1). 


Fig.  2.  (a)  A  schematic  diagram  showing  the  positions  of  fluorescence 
sphere  filled  with  ICG  in  the  Intralipid- 10%  suspension,  (b)  Fluorescence 
image  at  830-nm.  Independent  spatial  intensity  distributions 
(independent  components)  in  (c)  the  entrance,  and  (d)  exit  plane  of  the 
specimen. 


The  shape  of  the  4-mm  diameter  object 
was  reconstructed  by  using  back-projection 
Fourier  transforms.  Figure  3(c)  shows  the  pro¬ 
jection  onto  the  x-y  plane,  and  Figure  3(d) 
shows  the  corresponding  spatial  intensity  pro¬ 
files.  The  full  width  half  maximum  (FWHM)  of 
the  profile  in  the  x-direction  is  ~16-mm,  while 
that  in  the  y-direction  is  18  mm.  We  attribute 
this  change  in  shape  and  size  to  the  projection 
of  the  hook  being  superimposed  onto  that  of  the 
sphere.  This  artifact  was  much  reduced  for  the 
12-mm  diameter  sphere  where  the  shape  of  the 
hook  could  be  distinguished,  and  the  FWHM  of 
the  profiles  were  only  14  mm,  much  closer  to 
the  actual  value. 


5.  DISCUSSION 
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Fig.  3  (a)-(b)  Green’s  functions  (solid  line)  fits  to  Independent  spatial 
intensity  (circle)  from  Fig.  2  (c)-(d)  respectively,  (c)-(d)  Shape 
reconstructed  using  back-projection  Fourier  transform  with 
corresponding  spatial  intensity  profile 


The  approach  presented  in  this  article  effec¬ 
tively  used  a  multiple  source  excitation  (real¬ 
ized  through  scanning  of  the  sample  in  the  x-y 
plane)  and  multiple  detectors  (each  pixel  on 
the  CCD  may  be  viewed  as  a  detector)  for 
data  acquisition  providing  sufficient  data  for 
extracting  useful  location  information  about 
the  fluorescent  target.  While  OPTICA  can  be 
used  for  locating  absorptive  and  scattering 
targets  as  well,'^  the  advantage  for  the  fluo¬ 
rescent  target  is  that  the  background  light  can 
be  minimized  by  use  of  appropriate  filters, 
making  it  a  ‘zero-background’  measurement, 
as  opposed  to  the  situations  involving  absorp¬ 
tive  and  scattering  targets  when  changes  in 
high  light  levels  are  measured.  The  3-D  local¬ 
ization  of  the  object  is  achieved  with  a  high 
accuracy,  and  shorter  computation  time  com¬ 
pared  to  other  commonly  used  methods  for  3- 
D  inverse  image  reconstruction.  The  ability  to 
discern  the  ~4-mm  diameter  ICG  sphere  dem¬ 
onstrate  the  potential  to  detect  breast  tumors 
in  early  stage  of  growth  when  the  treatment  is 
most  effective. 
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ABSTRACT 

A  new  imaging  approach  for  three-dimensional  localization  and  characterization  of  absorptive,  scattering  or 
fluorescent  objects  in  a  turbid  medium  is  presented  and  demonstrated  using  simulated  and  experimental  data 
This  approach  uses  a  multi-source  and  multi-detector  signal  acquisition  scheme  and  independent  component 
analysis  (ICA)  from  information  theory  for  target  localization  and  characterization.  Independent  component 
analysis  of  the  perturbation  in  the  spatial  intensity  distribution  or  the  fluorescent  signal  measured  on  the  medium 
boundary  sorts  out  the  embedded  objects.  The  location  and  optical  characteristics  (size,  shape  and  optical 
property)  of  the  embedded  objects  are  obtained  from  a  Green’s  function  analysis  based  on  an  appropriate  model 
for  light  propagation  in  the  background  medium  and  back-projections  of  the  retrieved  independent  components. 

Keywords:  optical  imaging,  independent  component  analysis,  absorption,  scattering,  fluorescence 


1.  INTRODUCTION 

Optical  probing  of  the  interior  of  highly  scattering  colloidal  suspensions  and  biological  materials  has  attracted 
much  attention  over  the  last  decade.  In  particular,  biomedical  optical  tomography  and  spectroscopy  using  near 
infrared  (NIR)  light  has  been  actively  pursued  to  provide,  for  example,  functional  information  about  brain  ac¬ 
tivities  and  diagnostic  information  about  breast  tumor  because  NIR  light  can  penetrate  deep  into  tissues  and  is 
sensitive  to  certain  physiological  parameters  such  as  hemoglobin  concentration  and  blood  flow!^"^  Both  endoge¬ 
nous  contrasts  such  as  increased  hemoglobin  absorbance,  and  exogenous  contrasts  such  as  fluorescence  agents 
have  been  studied  for  the  purpose  of  medical  diagnostics  and  prognostic  measure. 

Most  works^’  in  imaging  of  endogenous  or  exogenous  contrasts  in  turbid  media  adopt  an  inverse  problem 
formalism  based  on  the  diffusion  approximation  or  radiative  transfer  of  light  propagation  in  highly  scattering 
media.  A  map  of  the  spatial  distribution  of  the  endogenous  contrast  inside  the  medium  is  reconstructed  by 
matching  the  detected  light  intensity  change  due  to  the  contrast  to  that  calculated  by  a  forward  model  for  light 
propagation  in  that  medium.  A  map  of  the  spatial  distribution  of  the  exogenous  contrast  is  reconstructed,  in  a 
similar  fashion,  by  matching  the  detected  fluorescence  to  that  calculated  by  a  forward  model  for  the  propagation 
of  the  fluorescence  light  excited  by  an  incident  light  beam  in  that  medium.  The  commonly  used  forward  models 
include  the  radiative  transfer  equation  (RTE),’^’  the  diffusion  approximation  (DA)  of  RTE,^  and  random  walk 
of  photons, 

Further  author  information:  (Send  correspondence  to  M.  X.) 

M.  X.;  Email:  minxu@sci.ccny.cuny.edu 


modi  i 


where  i  | 
is  the  c  I 
of  the  i  j 
spective  | 

-  I 


Optical  Tomography  and  Spectroscopy  of  Tissue  VI,  edited  by  Britton  Chance,  Robert  fl.  Alfafj^  j 
Bruce  J.  Tromberg,  Mamoru  Tamura,  Eva  M.  Sevick-Muraca,  Proc.  of  SPIE  Vol.  569a  j 
(SPIE,  Bellingham,  WA,2005)  •  1605-7422/05/$15  ■  doi:  10.1117/12.592153 


In  this  paper,  we  present  an  alternative  approach  using  independent  component  analysis  (ICA)  of  infor¬ 
mation  theory  to  detect  absorptive,  scattering,  or  fluorescent  inhomogeneities  within  flssue-Uke  turbid  media. 
We  use  a  multi-source  multi-detector  scheme  of  data  acquisition  where  Ught  intensities  are  measured  on  the 
boundary  of  tlie  medium  by  a  two  dimensional  detector  array  when  an  external  point  source  scans  on  the  other 
side  of  the  medium.  The  recorded  data  set  (after  removal  of  background  signal  in  imaging  of  absorption  or 
scattering  contrasts)  provides  a  set  of  observations  made  in  multiple  channels  (detectors)  from  virtual  sources 
(the  inhomogeneities  illuminated  by  the  incident  wave)  over  multiple  discrete  sampling  points  (the  source  scan¬ 
ning  positions).  Absorptive,  scattering  or  fluorescent  inhomogeneities  appear  as  statistically  independent  virtual 
sources.  These  virtual  sources  are  recovered  by  independent  component  analysis  of  the  data  set  and  sorts  out 
the  inhomogeneities.  The  location  and  characterization  (size,  shape  and  optical  property)  of  inhomogeneities 
are  obtained  from  a  Green’s  function  analysis  based  on  any  appropriate  model  of  light  propagation  in  the  back¬ 
ground  medium  and  back-projections  of  the  retrieved  independent  components.  This  approach  is  best  suited  for 
small  inhomogeneitieS. 
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We  refer  to  this  information  theory-inspired  approach  as  optical  imaging  using  independent  component 
analysis,  abbreviated  as,  OPTICA.  ICA  has  been  successfully  applied  in  various  other  applications  such  as  elec¬ 
troencephalogram  and  nuclear  magnetic  resonance  spectroscopy.^^"^®  The  novelty  of  OPTICA  over  other  ICA 
applications  is  that  OPTICA  associates  directly  the  independent  components  to  the  Green’s  functions  respon¬ 
sible  for  light  propagation  in  the  turbid  medium  from  the  inhomogeneities  to  the  source  and  the  detector,  and 
therefore  the  retrieved  independent  components  can  be  used  to  locate  and  characterize  the  inhomogeneities. 

2.  THEORY 

2.1.  Absorptive  and  scattering  inhomogeneities 

Absorptive  and  scattering  inhomogeneities  within  a  turbid  medium  alter  the  propagation  of  Ught  through  the 
medium.  The  alteration  of  the  spatial  distribution  of  the  Ught  intensity  on  the  boundary  of  the  medium  depends 
on  two  Green’s  functions  (propagators)  G(r,  r^,  lu)  and  G{T(1,  r,  ui)  describing  Ught  propagation  from  the  source 
at  to  the  inhomogeneity  at  r  and  from  the  inhomogeneity  to  the  detector  at  j^i,  respectively,  where  tu  is  the 
modulation  frequency  of  Ught, 

In  a  parallel  geometry,  the  scattered  wave  due  to  the  presence  of  J  absorptive  inhomogeneities  and  J 
scattering  inhoraogeneities  can  be  written  as’^® 


-^sca.{r  d,rs,uj)  =  '^G{rd,Tj,uj)qjG{rj,rs,uj)  +  Y^gzirj,rd,uj)q'jgz{Tj,rs,u) 


(1) 


j=i 


j=i 


J' 


+  M  ^dg±irj ,  rd,  u;)q'jPsj  cos  Osgxi^j,  r^,  w) 
J' 

+  (rj ,  Vd,  <^)q'jPsj  sin  OsPx (vj ,  ,  w) 


where  qj  =  5 p,a{Tj)cVj  is  the  absorption  strength  of  the  jth  inhomogeneity  of  volume  Vj,  q'j  =  5D{rj)cVj 
is  the  diffusion  strength  of  the  jth  scattering  inhomogeneity  of  volume  5pa  and  SD  are  the  deviation 
of  the  absorption  and  diffusion  coefficients  of  the  inhomogeneity  from  that  of  the  background  medium  re¬ 
spectively,  c  is  the  speed  of  Ught  in  the  background  medium,  pij  =  {xd  —  xj)^  -|-  {yd  —  VjY  and  Psj  — 

y  i^s  -  +  iVa  —  VjY  are  the  projections  of  -  Vj  and  —  rj  on  the  xy  plane,  respectively,  and  9d  and 
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Os  are  the  azimuthal  angles  of  -  Vj  and  -  Tj,  respectively.  The  expressions  of  the  Green’s  function  G  and 
the  two  auxiliary  functions  and  g_\_  depend  on  the  light  propagation  model  in  the  background  medium.^® 

In  OPTICA,  the  absorptive  or  scattering  inhomogeneities  illuminated  by  the  incident  wave  are  assumed  to  be 
virtual  sources,  and  the  perturbation  of  the  spatial  distribution  of  the  light  intensity  on  the  boundary  described  by 
Eq.  (1)  may  be  considered  to  be  a  linear  mixture  of  signals  arriving  from  these  virtual  sources.  One  absorptive 
inhomogeneity  is  represented  by  one  virtual  source  rs)  with  a  mixing  vector  Gira,  rj).  One  scattering 

inhomogeneity  is  represented  by  three  virtual  sources: 


with  mixing  vectors 


qjg^(Tj,rs),  q'jPsj  CQsOsgi.{rj,Ts),  (/jpsj  sin0s5x(rj, r^), 
gz{Tj,Td),  prfjCos^d9x(rj,rd),  PdjsmOdgi.{Tj,rd) 


respectively. 

Owing  to  the  statistical  independence  between  virtual  sources,  independent  component  analysis  of  the 
recorded  data  set  can  identify  prominent  independent  components.  Each  independent  component  consists  of 
one  virtual  source  and  its  corresponding  weighting  vector.  The  number  of  leading  independent  components  is 
same  as  the  number  of  inhomogeneities.  The  details  of  the  theoretical  formalism  of  OPTICA  for  absorption  or 
scattering  inhomogeneities  were  presented  elsewhere.^® 

2.2.  Fluorescent  inhomogeneities 

Light  propagation  in  a  highly  scattering  medium  with  fluorescent  targets  excited  by  an  external  source  is  approxi¬ 
mately  described  by  coupled  diffusion  equations  at  the  excitation  and  emission  wavelengths?’  ®  The  fluorescence 
signal  Um{rd,  r®,  w)  can  be  expressed  in  terms  of  the  two  Green’s  functions  Gx{t,  Tg,  uj)  and  Gmi^d,  u>)  de¬ 
scribing  light  propagation  from  the  source  at  to  the  target  at  r  at  the  excitation  wavelength  and  the  light 
propagation  from  the  target  to  the  detector  at  at  the  emission  wavelength  A^.  The  subscripts  x  and  m  refer 
to  the  quantities  associated  with  the  excitation  and  emission  wavelengths,  respectively. 

Assuming  the  jth  fluorophores  is  contained  in  volume  Vj  centered  at  rj,  the  fluorescence  signal  under 
illumination  by  a  unit  point  source  at  Tg  is  given  by 

j 

j=t 

where  fj{uj)  =  7(rj)cl^/  [1  —  iLJT{Tj)]  represents  the  fluorescence  strength  of  the  jth  inhomogeneity,  7  is  the 
fluorescent  yield,  and  r  is  the  fluorescence  lifetime. 

Eq.  (4)  takes  exactly  the  same  form  of  the  scattered  wave  Eq.  (1)  for  absorptive  inhomogeneities  only  (set¬ 
ting  g'  =  0inEq.(l)).  The  yth  virtual  source  ^Ga;(rj,rs,  w)  is  weighted  by  a  mixing  vector  G^(r(i,  rj-.w)  and 
all  the  weighted  virtual  sources  are  mixed  to  produce  the  detected  signal  Uni^d,  Ts,ui).  By  seeking  the  maxi¬ 
mal  mutual  independence,  both  virtual  sources  and  mixing  vectors  can  be  obtained  by  independent  component 
analysis  of  the  multi-source  multi-detector  data  set  Um{ttd>  fs)  w)- 

Let  us  denote  the  recovered  jth  virtual  source  by  %(rg)  oc  Gx{fj,  rg,a;)  and  the  recovered  corresponding 
mixing  vector  by  oc  Gm{ttd:  ry,  uj).  Both  localization  and  strength  of  the  jth  object  can  be  computed  by 

a  least  square  fitting  procedure; 


-  Gxirj,  Tg)]  +  Y,  \Pj  ^oy(rd)  -  Gm{rd,  ry)l 

I  r.  ;  :  rj  ^ 


II 

sponc 
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The  fitting  yields  the  location  r^-  of  and  the  two  scaling  constants  aj  and  Pj  for  the  jth  inhomogeneity.  The 
fluorescence  strength  is  then 

The  size  and  shape  of  the  ii^omogeheity  can  be  estimated  after  the  fluorescent  object  has  been  locahzed. 
The  fluorescence  signal  due  &lh©'jth  object  is  (r(i,  r^,  w)  =  aj  (rrf)5j  (r^).  One  can  estimate  the  volume 
Vj  of  the  jth  object  by  its  thickness  profile  A.^(p)  in  the  z  direction  where  p  is  the  lateral  coordinate  on  the 
z  =  Zj  plane  with  zy  the  «  coordinate  of  the  center  of  die  yth  object.  The  fluorescence  signal  from  the  jith  object 
is  approximately  given  by: 

Umj{rd, J Gm{pd  -  P,  tu)Azj(p)6V(p  -  Ps, i^]dp  (6) 

where  the  integration  is  performed  over  the  z  =  Zj  plane,  and  pd  and  p,,  are  the  lateral  coordinates  of  the 
detector  and  the  source  respectively.  In  the  Fourier  space,  the  thickness  A.5  crui  be  solved  from  Eq.  (6)  and  is 
given  by 


A^ry(q) 


/jH  Gm(q  -  qs,a;)G*(qs,  w) ' 


(7) 


where  q  and  q^  are  the  spatial  frequency  on  the  xy  plane  and*  denotes  complex  conjugate.  The  inverse  Fourier 
transform  of  Azyfq)  yields  the  thickness  profile  of  the  fluorescent  target  in  the  2  direction.  The  FWHM  (full 
width  at  half  maximum  value)  and  the  contour  of  the  thickness  profile  provide  an  estimation  of  size  and  shape 
of  the  j'th  target,  respectively.  This  procedure  to  obtain  size  and  shape  of  the  fluorescent  inhomogeneity  also 
applies  to  the  absorptive  and  scattering  inhomogeneities. 


3.  RESULTS 

The  first  test  of  OPTICA  involved  imaging  of  two  absorptive  inhomogeneities  using  simulated  data.  In  simu¬ 
lation,  the  absorptive  objects  are  placed  at  (50, 60, 20)mm  and  (30, 30, 30)mm,  respectively,  inside  a  slab  of 
volume  90x90x50mm2  shown  schematically  in  Fig.  1(a).  Each  inhomogeneity  has  a  volume  of  6.75  x  6.75  x 
4.69mm^  and  an  absorption  strength  of  unity.  The  absorption  and  diffusion  coefficients  of  the  uniform  slab  is 
Pa  =  l/300imn“^and  D  =  l/3mm  respectively,  close  to  that  of  human  breast  tissue.^'^  The  incident  CWbeam 
scans  a  set  of  21  x  21  grid  points  covering  an  area  of  90  x  90mrE?.  The  spacing  between  two  consecutive  grid 
points  is  4.5mm.  The  light  intensity  on  the  other  side  of  the  slab  is  recorded  by  a  CCD  camera  on  42  x  42  grid 
points  covering  the  same  area.  10%  Gaussian  noise  was  added  to  the  simulated  light  intensity  change  on  the 
detector  plane. 

Independent  component  analysis  of  the  perturbation  in  the  spatial  intensity  distributions  provided  corre¬ 
sponding  independent  intensity  distributions  on  the  source  and  detector  planes.  IC A  generated  independent 
intensity  distributions  on  the  source  and  detector  planes  are  shown  in  Fig.  1(b),  for  the  two  absorptive  inhomo¬ 
geneities.  Locations  of  the  absorptive  objects  are  obtained  from  fitting  these  independent  component  intensity 
distributions  to  those  of  the  diffusion  approximation  in  a  slab  by  least  square.  The  first  object  is  found  at 
(50.0, 60.0, 21.1)mm  and  the  second  one  at  (30.0, 30.0, 29.9)mm.  The  xy  coordinates  of  both  objects  agree 
exactly  with  their  known  values  and  the  z  coordinate  of  both  objects  is  within  1.2mm  of  their  known  center 
locations.  The  strengths  of  the  two  objects  are  ql  =  0.88  and  q2  =  0.88  respectively,  with  an  error  of  ~  12%  of 
the  true  values. 

The  thickness  map  for  the  two  absorptive  objects  is  back-projected  using  a  formula  similar  to  Eq.  (7)  and 
presented  in  Fig.  2.  For  each  case,  the  Green’s  function  at  the  source  and  detector  planes  (circle)  are  shown  with 
the  fitting  in  solid  fines  on  the  first  row.  The  longitudinal  thickness  map  of  the  target  and  the  thickness  profiles 
along  A  and  Y  directions  are  also  displayed  on  the  second  row.  Both  targets  are  observed  to  have  approximately 
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Figure  1.  (a)  The  light  intensity  on  one  side  of  the  slab  is  measured  when  a  point  source  scans  on  the  other  side.  Two 
inhomogeneities  are  placed  at  (50,60, 20)mm  and  (30, 30,30)mm  inside  a  90  x  90  x  50mm®  slab,  (b)  Normalized 
independent  spatial  intensity  distributions  on  the  source  plane  (the  first  row),  the  detector  plane  (the  second  row),  and  the 
least  square  fitting  (the  third  row).  The  left  column  is  for  the  first  absorptive  inhomogeneity  at  (50, 60, 20)mm  and  the 
right  column  is  for  the  second  absorptive  inhomogeneity  at  (30, 30, 30)mm.  On  the  third  row,  the  horizontal  profile  of 
intensity  distributions  on  the  source  plane  (diamond)  and  on  the  detector  plane  (circle)  are  displayed,  and  solid  lines  show 
the  respective  Green’s  function  fit  used  for  obtaining  locations  and  strengths  of  objects.  The  noise  level  is  10%, 

a  square  shape.  The  FWHM  of  the  peak  for  the  first  and  the  second  target  are  found  to  be  16.0mm  and  13.6mm, 
and  9.5mm  and  10.5mm  in  the  x  and  y  directions,  respectively.  These  values  should  be  compared  to  the  size  of 
the  absorptive  target  of  6.75mm  in  the  xy  plane.  The  estimation  of  the  size  and  shape  of  the  second  absorptive 
object  is  more  accurate  as  it  is  closer  to  the  detector  plane  where  the  grid  spacing  is  finer. 


The  second  example  is  to  image  a  spherical  fluorescent  target  placed  inside  a  250mmx250mm  x  50ram 
slab  filled  with  Intralipid- 10%  aqueous  suspension  using  experimental  data  [see  Fig.  3(a)].  The  fluorescent  tar¬ 
get  was  a  9.0mm  diameter  sphere  filled  with  a  solution  of  Indocyanine  green  (ICG)  dye  in  water  that  could  be 
excited  in  the  650nm-800nm  spectral  range.  The  wavelength  of  the  CW  incident  light  is  =  785nm.  Two  long 

wavelength  pass  absorption  filters  were  placed  before  the  CCD  camera  to  block  the  excitation  wavelength  and 
allow  the  fluorescence  light  to  pass.  The  wavelength  of  the  peak  fluorescence  light  adjusted  by  the  filtering  and 
the  camera  response  efficiency  is  about  \n  =  870nm.  The  Intralipid- 10%  suspension  is  diluted  with  pure  water 
such  that  the  transport  mean  free  paths  and  absorption  coefficients  are  =  1.01mm  and  fiax  —  0.0022mm 
at  the  excitation  wavelength  and  =  1.14mm  and  fj-am  ==  0.0054mm~^  at  the  emission  wavelength,  respec¬ 
tively.  The  sphere  is  placed  at  the  position  (11, 9, 32)mm  inside  the  slab.  The  sphere  is  33mm  away  from  the 
input  window.  The  point  source  scans  over  a  10  x  10  grid  system  with  spacing  2.5mm  between  consecutive 
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(a)  (b) 

Figure  2,  Fitting  of  the  independent  intensity  distribution  of  the  absorptive  objects  of  volume  6.75  x  6.75  x  4.69mm  ^  inside 
the  slab  described  in  the  text  for  (a)  the  first  absorptive  object  at  (50, 60, 20)mm  and  (b)  the  second  absorptive  object  at 
(30, 30, 30)mm.  For  each  case,  the  Green’s  function  at  the  source  and  detector  planes  (circle)  are  shown  with  the  fitting  in 
solid  lines  [the  first  row].  The  longitudinal  thickness  map  of  the  target  and  the  thickness  profiles  along  X  and  Y  directions 
are  also  displayed  [the  second  row]. 


I  10 


40t  - 
-20 


5  10  15  20  25 

:  X  (mmV  : 


(b) 

Figure  3.  (a)  The  schematic  diagram  showing  a  fluorescent  sphere  embedded  inside  a  slab,  (b)  Independent  intensity 
distributions  of  the  fluorescence  from  the  target  generated  by  ICA  at  the  detector  plane  (Left)  and  at  the  source  plane 
(Right).  The  fluorescent  sphere  is  located  on  the  2r  =  32mm  plane. 

giidsv  The  grid  spacing  on  the  detector  plane  is  1.23mm. 

Independent  intensity  distributions  at  the  detector  plane  and  the  source  plane  obtained  by  ICA  is  shown  in 
Fig.  3(b).  The  fluorescent  target  is  found  to  be  centered  at  {x,y)  =  (ll,9)mm  and  2:  =  33mm  away  from 
the  input  window  by  fitting  independent  intensity  distributions  at  the  detector  plane  and  the  source  plane  to  the 
Green’s  functions  in  the  diffusion  approximation  [see  Fig.  4(a)  and  (b)].  The  (x,y)  coordinate  agrees  exactly 
with  the  known  xy  coordinates  and  the  z  coordinate  is  within  1mm  error  of  the  known  z  coordinate.  The 
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Figure  4.  Fitting  of  the  independent  intensity  distribution  of  fluorescence  from  a  sphere  of  diameter  9inm  embedded  in 
Intralipid- 10%  solution  to  the  Green’s  function  (a)  at  the  detector  plane  and  (b)  at  the  source  plane,  (c)  The  longitudinal 
thickness  map  of  the  target  centered  at  (11, 9)mm.  (d)  The  thickness  profiles  along  X  and  Y  directions. 

fluorescence  strength  is  obtained  to  be  /  =  0.0462mir]? /ps. 

The  thickness  map  is  back-projected  using  Eq.  (7)  and  presented  in  Fig.  4(c).  The  horizontal  and  vertical 
thickness  profiles  of  A^/A^max  passing  through  the  center  are  also  plotted  in  Fig.  4(d).  The  target  is  observed 
to  have  a  circular  shape.  The  FWHM  of  the  peak  found  to  be  d  -  11.5mm.  This  value  should  be  compared  to 
the  diameter  of  the  fluorescent  target  9mm. 
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